c763e2e762b0ea2cdb364ea26fd767d1ad6be6e0
[occt.git] / src / BRepLib / BRepLib.cxx
1 // Created on: 1993-12-15
2 // Created by: Remi LEQUETTE
3 // Copyright (c) 1993-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 //pmn 26/09/97 Add parameters of approximation in BuildCurve3d
18 //  Modified by skv - Thu Jun  3 12:39:19 2004 OCC5898
19
20 #include <Adaptor3d_CurveOnSurface.hxx>
21 #include <AdvApprox_ApproxAFunction.hxx>
22 #include <AppParCurves_MultiBSpCurve.hxx>
23 #include <AppParCurves_MultiCurve.hxx>
24 #include <Approx_CurvilinearParameter.hxx>
25 #include <Approx_SameParameter.hxx>
26 #include <Bnd_Box.hxx>
27 #include <BRep_Builder.hxx>
28 #include <BRep_CurveRepresentation.hxx>
29 #include <BRep_GCurve.hxx>
30 #include <BRep_ListIteratorOfListOfCurveRepresentation.hxx>
31 #include <BRep_ListOfCurveRepresentation.hxx>
32 #include <BRep_TEdge.hxx>
33 #include <BRep_TFace.hxx>
34 #include <BRep_Tool.hxx>
35 #include <BRep_TVertex.hxx>
36 #include <BRepAdaptor_HCurve.hxx>
37 #include <BRepAdaptor_HCurve2d.hxx>
38 #include <BRepAdaptor_HSurface.hxx>
39 #include <BRepAdaptor_Surface.hxx>
40 #include <BRepBndLib.hxx>
41 #include <BRepClass3d_SolidClassifier.hxx>
42 #include <BRepLib.hxx>
43 #include <BSplCLib.hxx>
44 #include <ElSLib.hxx>
45 #include <Extrema_LocateExtPC.hxx>
46 #include <GCPnts_QuasiUniformDeflection.hxx>
47 #include <Geom2d_BSplineCurve.hxx>
48 #include <Geom2d_Curve.hxx>
49 #include <Geom2d_TrimmedCurve.hxx>
50 #include <Geom2dAdaptor.hxx>
51 #include <Geom2dAdaptor_Curve.hxx>
52 #include <Geom2dAdaptor_HCurve.hxx>
53 #include <Geom2dConvert.hxx>
54 #include <Geom_BSplineCurve.hxx>
55 #include <Geom_BSplineSurface.hxx>
56 #include <Geom_Curve.hxx>
57 #include <Geom_Plane.hxx>
58 #include <Geom_RectangularTrimmedSurface.hxx>
59 #include <Geom_Surface.hxx>
60 #include <Geom_TrimmedCurve.hxx>
61 #include <GeomAdaptor_Curve.hxx>
62 #include <GeomAdaptor_HCurve.hxx>
63 #include <GeomAdaptor_HSurface.hxx>
64 #include <GeomAdaptor_Surface.hxx>
65 #include <GeomLib.hxx>
66 #include <GeomLProp_SLProps.hxx>
67 #include <gp.hxx>
68 #include <gp_Ax2.hxx>
69 #include <gp_Pln.hxx>
70 #include <Poly_PolygonOnTriangulation.hxx>
71 #include <Poly_Triangulation.hxx>
72 #include <Precision.hxx>
73 #include <ProjLib_ProjectedCurve.hxx>
74 #include <Standard_ErrorHandler.hxx>
75 #include <Standard_Real.hxx>
76 #include <TColgp_Array1OfPnt.hxx>
77 #include <TColgp_Array1OfPnt2d.hxx>
78 #include <TColStd_Array1OfReal.hxx>
79 #include <TColStd_MapOfTransient.hxx>
80 #include <TopExp.hxx>
81 #include <TopExp_Explorer.hxx>
82 #include <TopoDS.hxx>
83 #include <TopoDS_Edge.hxx>
84 #include <TopoDS_Face.hxx>
85 #include <TopoDS_Shape.hxx>
86 #include <TopoDS_Solid.hxx>
87 #include <TopoDS_Vertex.hxx>
88 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
89 #include <TopTools_ListIteratorOfListOfShape.hxx>
90 #include <TopTools_MapOfShape.hxx>
91 #include <TShort_HArray1OfShortReal.hxx>
92 #include <TColgp_Array1OfXY.hxx>
93
94 #include <algorithm>
95
96 // TODO - not thread-safe static variables
97 static Standard_Real thePrecision = Precision::Confusion();     
98 static Handle(Geom_Plane) thePlane;
99
100 //=======================================================================
101 // function: BRepLib_ComparePoints
102 // purpose:  implementation of IsLess() function for two points
103 //=======================================================================
104 struct BRepLib_ComparePoints {
105   bool operator()(const gp_Pnt& theP1, const gp_Pnt& theP2)
106   {
107     for (Standard_Integer i = 1; i <= 3; ++i) {
108       if (theP1.Coord(i) < theP2.Coord(i)) {
109         return Standard_True;
110       }
111       else if (theP1.Coord(i) > theP2.Coord(i)) {
112         return Standard_False;
113       }
114     }
115     return Standard_False;
116   }
117 };
118
119
120 //=======================================================================
121 //function : Precision
122 //purpose  : 
123 //=======================================================================
124
125 void BRepLib::Precision(const Standard_Real P)
126 {
127   thePrecision = P;
128 }
129
130 //=======================================================================
131 //function : Precision
132 //purpose  : 
133 //=======================================================================
134
135 Standard_Real  BRepLib::Precision()
136 {
137   return thePrecision;
138 }
139
140 //=======================================================================
141 //function : Plane
142 //purpose  : 
143 //=======================================================================
144
145 void  BRepLib::Plane(const Handle(Geom_Plane)& P)
146 {
147   thePlane = P;
148 }
149
150
151 //=======================================================================
152 //function : Plane
153 //purpose  : 
154 //=======================================================================
155
156 const Handle(Geom_Plane)&  BRepLib::Plane()
157 {
158   if (thePlane.IsNull()) thePlane = new Geom_Plane(gp::XOY());
159   return thePlane;
160 }
161 //=======================================================================
162 //function : CheckSameRange
163 //purpose  : 
164 //=======================================================================
165
166 Standard_Boolean  BRepLib::CheckSameRange(const TopoDS_Edge& AnEdge,
167   const Standard_Real Tolerance) 
168 {
169   Standard_Boolean  IsSameRange = Standard_True,
170     first_time_in = Standard_True ;
171
172   BRep_ListIteratorOfListOfCurveRepresentation an_Iterator
173     ((*((Handle(BRep_TEdge)*)&AnEdge.TShape()))->ChangeCurves());
174
175   Standard_Real first, last;
176   Standard_Real current_first =0., current_last =0. ;
177   Handle(BRep_GCurve) geometric_representation_ptr ;
178
179   while (IsSameRange && an_Iterator.More()) {
180     geometric_representation_ptr =
181       Handle(BRep_GCurve)::DownCast(an_Iterator.Value());
182     if (!geometric_representation_ptr.IsNull()) {
183
184       first = geometric_representation_ptr->First();
185       last =  geometric_representation_ptr->Last();
186       if (first_time_in ) {
187         current_first = first ;
188         current_last = last   ;
189         first_time_in = Standard_False ;
190       }
191       else {
192         IsSameRange = (Abs(current_first - first) <= Tolerance) 
193           && (Abs(current_last -last) <= Tolerance ) ;
194       }
195     }
196     an_Iterator.Next() ;
197   }
198   return IsSameRange ;
199 }
200
201 //=======================================================================
202 //function : SameRange
203 //purpose  : 
204 //=======================================================================
205
206 void BRepLib::SameRange(const TopoDS_Edge& AnEdge,
207   const Standard_Real Tolerance) 
208 {
209   BRep_ListIteratorOfListOfCurveRepresentation an_Iterator
210     ((*((Handle(BRep_TEdge)*)&AnEdge.TShape()))->ChangeCurves());
211
212   Handle(Geom2d_Curve) Curve2dPtr, Curve2dPtr2, NewCurve2dPtr, NewCurve2dPtr2;
213   TopLoc_Location LocalLoc ;
214
215   Standard_Boolean first_time_in = Standard_True,
216     has_curve,
217     has_closed_curve ;
218   Handle(BRep_GCurve) geometric_representation_ptr ;
219   Standard_Real first,
220     current_first,
221     last,
222     current_last ;
223
224   const Handle(Geom_Curve) C = BRep_Tool::Curve(AnEdge,
225     LocalLoc,
226     current_first,
227     current_last);
228   if (!C.IsNull()) {
229     first_time_in = Standard_False ;
230   }
231
232   while (an_Iterator.More()) {
233     geometric_representation_ptr =
234       Handle(BRep_GCurve)::DownCast(an_Iterator.Value());
235     if (! geometric_representation_ptr.IsNull()) {
236       has_closed_curve =
237         has_curve = Standard_False ;
238       first = geometric_representation_ptr->First();
239       last =  geometric_representation_ptr->Last();
240       if (geometric_representation_ptr->IsCurveOnSurface()) {
241         Curve2dPtr = geometric_representation_ptr->PCurve() ; 
242         has_curve = Standard_True ;
243       }
244       if (geometric_representation_ptr->IsCurveOnClosedSurface()) {
245         Curve2dPtr2 = geometric_representation_ptr->PCurve2() ;
246         has_closed_curve = Standard_True ;
247       }
248       if (has_curve || has_closed_curve) {
249         if (first_time_in) {
250           current_first = first ;
251           current_last = last ;
252           first_time_in = Standard_False ;
253         }
254
255         if (Abs(first - current_first) > Precision::Confusion() ||
256           Abs(last - current_last) > Precision::Confusion() )
257         {
258           if (has_curve)
259           {
260             GeomLib::SameRange(Tolerance,
261               Curve2dPtr,
262               geometric_representation_ptr->First(),
263               geometric_representation_ptr->Last(),
264               current_first,
265               current_last,
266               NewCurve2dPtr);
267             geometric_representation_ptr->PCurve(NewCurve2dPtr) ;
268           }
269           if (has_closed_curve)
270           {
271             GeomLib::SameRange(Tolerance,
272               Curve2dPtr2,
273               geometric_representation_ptr->First(),
274               geometric_representation_ptr->Last(),
275               current_first,
276               current_last,
277               NewCurve2dPtr2);
278             geometric_representation_ptr->PCurve2(NewCurve2dPtr2) ;
279           }
280         }
281       }
282     }
283     an_Iterator.Next() ;
284   }
285   BRep_Builder B;
286   B.Range(TopoDS::Edge(AnEdge),
287     current_first,
288     current_last) ;
289
290   B.SameRange(AnEdge,
291     Standard_True) ;
292 }
293
294 //=======================================================================
295 //function : EvaluateMaxSegment
296 //purpose  : return MaxSegment to pass in approximation, if MaxSegment==0 provided
297 //=======================================================================
298
299 static Standard_Integer evaluateMaxSegment(const Standard_Integer aMaxSegment,
300   const Adaptor3d_CurveOnSurface& aCurveOnSurface)
301 {
302   if (aMaxSegment != 0) return aMaxSegment;
303
304   Handle(Adaptor3d_HSurface) aSurf   = aCurveOnSurface.GetSurface();
305   Handle(Adaptor2d_HCurve2d) aCurv2d = aCurveOnSurface.GetCurve();
306
307   Standard_Real aNbSKnots = 0, aNbC2dKnots = 0;
308
309   if (aSurf->GetType() == GeomAbs_BSplineSurface) {
310     Handle(Geom_BSplineSurface) aBSpline = aSurf->BSpline();
311     aNbSKnots = Max(aBSpline->NbUKnots(), aBSpline->NbVKnots());
312   }
313   if (aCurv2d->GetType() == GeomAbs_BSplineCurve) {
314     aNbC2dKnots = aCurv2d->NbKnots();
315   }
316   Standard_Integer aReturn = (Standard_Integer) (  30 + Max(aNbSKnots, aNbC2dKnots) ) ;
317   return aReturn;
318 }
319
320 //=======================================================================
321 //function : BuildCurve3d
322 //purpose  : 
323 //=======================================================================
324
325 Standard_Boolean  BRepLib::BuildCurve3d(const TopoDS_Edge& AnEdge,
326   const Standard_Real Tolerance,
327   const GeomAbs_Shape Continuity,
328   const Standard_Integer MaxDegree,
329   const Standard_Integer MaxSegment)
330 {
331   Standard_Integer //ErrorCode,
332     //                   ReturnCode = 0,
333     ii,
334     //                   num_knots,
335     jj;
336
337   TopLoc_Location LocalLoc,L[2],LC;
338   Standard_Real f,l,fc,lc, first[2], last[2],
339     tolerance,
340     max_deviation,
341     average_deviation ;
342   Handle(Geom2d_Curve) Curve2dPtr, Curve2dArray[2]  ;
343   Handle(Geom_Surface) SurfacePtr, SurfaceArray[2]  ;
344
345   Standard_Integer not_done ;
346   // if the edge has a 3d curve returns true
347
348
349   const Handle(Geom_Curve) C = BRep_Tool::Curve(AnEdge,LocalLoc,f,l);
350   if (!C.IsNull()) 
351     return Standard_True;
352   //
353   // this should not exists but UpdateEdge makes funny things 
354   // if the edge is not same range 
355   //
356   if (! CheckSameRange(AnEdge,
357     Precision::Confusion())) {
358       SameRange(AnEdge,
359         Tolerance) ;
360   }
361
362
363
364   // search a curve on a plane
365   Handle(Geom_Surface) S;
366   Handle(Geom2d_Curve) PC;
367   Standard_Integer i = 0;
368   Handle(Geom_Plane) P;
369   not_done = 1 ;
370
371   while (not_done) {
372     i++;
373     BRep_Tool::CurveOnSurface(AnEdge,PC,S,LocalLoc,f,l,i);
374     Handle(Geom_RectangularTrimmedSurface) RT = 
375       Handle(Geom_RectangularTrimmedSurface)::DownCast(S);
376     if ( RT.IsNull()) {
377       P = Handle(Geom_Plane)::DownCast(S);
378     }
379     else {
380       P = Handle(Geom_Plane)::DownCast(RT->BasisSurface());
381     }
382     not_done = ! (S.IsNull() || !P.IsNull()) ;
383   }
384   if (! P.IsNull()) {
385     // compute the 3d curve
386     gp_Ax2 axes = P->Position().Ax2();
387     Handle(Geom_Curve) C3d = GeomLib::To3d(axes,PC);
388     if (C3d.IsNull())
389       return Standard_False;
390     // update the edge
391     Standard_Real First, Last;
392
393     BRep_Builder B;
394     B.UpdateEdge(AnEdge,C3d,LocalLoc,0.0e0);
395     BRep_Tool::Range(AnEdge, S, LC, First, Last);
396     B.Range(AnEdge, First, Last); //Do not forget 3D range.(PRO6412)
397   }
398   else {
399     //
400     // compute the 3d curve using existing surface
401     //
402     fc = f ;
403     lc = l ;
404     if (!BRep_Tool::Degenerated(AnEdge)) {
405       jj = 0 ;
406       for (ii = 0 ; ii < 3 ; ii++ ) {
407         BRep_Tool::CurveOnSurface(TopoDS::Edge(AnEdge),
408           Curve2dPtr,
409           SurfacePtr,
410           LocalLoc,
411           fc,
412           lc,
413           ii) ;
414
415         if (!Curve2dPtr.IsNull() && jj < 2){
416           Curve2dArray[jj] = Curve2dPtr ;
417           SurfaceArray[jj] = SurfacePtr ;
418           L[jj] = LocalLoc ;
419           first[jj] = fc ;
420           last[jj] = lc ;
421           jj += 1 ;
422         }
423       }
424       f = first[0] ;
425       l = last[0] ;
426       Curve2dPtr = Curve2dArray[0] ;
427       SurfacePtr = SurfaceArray[0] ;
428
429       Geom2dAdaptor_Curve     AnAdaptor3dCurve2d (Curve2dPtr, f, l) ;
430       GeomAdaptor_Surface     AnAdaptor3dSurface (SurfacePtr) ;
431       Handle(Geom2dAdaptor_HCurve) AnAdaptor3dCurve2dPtr =
432         new Geom2dAdaptor_HCurve(AnAdaptor3dCurve2d) ;
433       Handle(GeomAdaptor_HSurface) AnAdaptor3dSurfacePtr =
434         new GeomAdaptor_HSurface (AnAdaptor3dSurface) ;
435       Adaptor3d_CurveOnSurface  CurveOnSurface( AnAdaptor3dCurve2dPtr,
436         AnAdaptor3dSurfacePtr) ;
437
438       Handle(Geom_Curve) NewCurvePtr ;
439
440       GeomLib::BuildCurve3d(Tolerance,
441         CurveOnSurface,
442         f,
443         l,
444         NewCurvePtr,
445         max_deviation,
446         average_deviation,
447         Continuity,
448         MaxDegree,
449         evaluateMaxSegment(MaxSegment,CurveOnSurface)) ;
450       BRep_Builder B;   
451       tolerance = BRep_Tool::Tolerance(AnEdge) ;
452       //Patch
453       //max_deviation = Max(tolerance, max_deviation) ;
454       max_deviation = Max( tolerance, Tolerance );
455       if (NewCurvePtr.IsNull())
456         return Standard_False;
457       B.UpdateEdge(TopoDS::Edge(AnEdge),
458         NewCurvePtr,
459         L[0],
460         max_deviation) ;
461       if (jj == 1 ) {
462         //
463         // if there is only one curve on surface attached to the edge
464         // than it can be qualified sameparameter
465         //
466         B.SameParameter(TopoDS::Edge(AnEdge),
467           Standard_True) ;
468       }
469     }
470     else {
471       return Standard_False ;
472     }
473
474   }         
475   return Standard_True;
476 }
477 //=======================================================================
478 //function : BuildCurves3d
479 //purpose  : 
480 //=======================================================================
481
482 Standard_Boolean  BRepLib::BuildCurves3d(const TopoDS_Shape& S) 
483
484 {
485   return BRepLib::BuildCurves3d(S,
486     1.0e-5) ;
487 }
488
489 //=======================================================================
490 //function : BuildCurves3d
491 //purpose  : 
492 //=======================================================================
493
494 Standard_Boolean  BRepLib::BuildCurves3d(const TopoDS_Shape& S,
495   const Standard_Real Tolerance,
496   const GeomAbs_Shape Continuity,
497   const Standard_Integer MaxDegree,
498   const Standard_Integer MaxSegment)
499 {
500   Standard_Boolean boolean_value,
501     ok = Standard_True;
502   TopTools_MapOfShape a_counter ;
503   TopExp_Explorer ex(S,TopAbs_EDGE);
504
505   while (ex.More()) {
506     if (a_counter.Add(ex.Current())) {
507       boolean_value = 
508         BuildCurve3d(TopoDS::Edge(ex.Current()),
509         Tolerance, Continuity,
510         MaxDegree, MaxSegment);
511       ok = ok && boolean_value ;
512     }
513     ex.Next();
514   }
515   return ok;
516 }
517 //=======================================================================
518 //function : UpdateEdgeTolerance
519 //purpose  : 
520 //=======================================================================
521
522 Standard_Boolean  BRepLib::UpdateEdgeTol(const TopoDS_Edge& AnEdge,
523   const Standard_Real MinToleranceRequested,
524   const Standard_Real MaxToleranceToCheck)
525 {
526
527   Standard_Integer curve_on_surface_index,
528     curve_index,
529     not_done,
530     has_closed_curve,
531     has_curve,
532     jj,
533     ii,
534     geom_reference_curve_flag = 0,
535     max_sampling_points = 90,
536     min_sampling_points = 30 ;
537
538   Standard_Real factor = 100.0e0,
539     //     sampling_array[2],
540     safe_factor = 1.4e0,
541     current_last,
542     current_first,
543     max_distance,
544     coded_edge_tolerance,
545     edge_tolerance = 0.0e0 ;
546   Handle(TColStd_HArray1OfReal) parameters_ptr ;
547   Handle(BRep_GCurve) geometric_representation_ptr ;
548
549   if (BRep_Tool::Degenerated(AnEdge)) return Standard_False ;
550   coded_edge_tolerance = BRep_Tool::Tolerance(AnEdge) ;
551   if (coded_edge_tolerance > MaxToleranceToCheck) return Standard_False ;
552
553   const Handle(BRep_TEdge)& TE = *((Handle(BRep_TEdge)*)&AnEdge.TShape());
554   BRep_ListOfCurveRepresentation& list_curve_rep = TE->ChangeCurves() ;
555   BRep_ListIteratorOfListOfCurveRepresentation an_iterator(list_curve_rep),
556     second_iterator(list_curve_rep) ;
557   Handle(Geom2d_Curve) curve2d_ptr, new_curve2d_ptr;
558   Handle(Geom_Surface) surface_ptr ;
559   TopLoc_Location local_location ;
560   GCPnts_QuasiUniformDeflection  a_sampler ;
561   GeomAdaptor_Curve  geom_reference_curve ;
562   Adaptor3d_CurveOnSurface  curve_on_surface_reference ; 
563   Handle(Geom_Curve) C = BRep_Tool::Curve(AnEdge,
564     local_location,
565     current_first,
566     current_last);
567   curve_on_surface_index = -1 ;
568   if (!C.IsNull()) {
569     if (! local_location.IsIdentity()) {
570       C = Handle(Geom_Curve)::
571         DownCast(C-> Transformed(local_location.Transformation()) ) ;
572     }
573     geom_reference_curve.Load(C) ;
574     geom_reference_curve_flag = 1 ;
575     a_sampler.Initialize(geom_reference_curve,
576       MinToleranceRequested * factor,
577       current_first,
578       current_last) ;
579   }
580   else {
581     not_done = 1 ;
582     curve_on_surface_index = 0 ;  
583
584     while (not_done && an_iterator.More()) {
585       geometric_representation_ptr =
586         Handle(BRep_GCurve)::DownCast(second_iterator.Value());
587       if (!geometric_representation_ptr.IsNull() 
588         && geometric_representation_ptr->IsCurveOnSurface()) {
589           curve2d_ptr = geometric_representation_ptr->PCurve() ;
590           local_location = geometric_representation_ptr->Location() ;
591           current_first = geometric_representation_ptr->First();
592           //first = geometric_representation_ptr->First();
593           current_last =  geometric_representation_ptr->Last();
594           // must be inverted 
595           //
596           if (! local_location.IsIdentity() ) {
597             surface_ptr = Handle(Geom_Surface)::
598               DownCast( geometric_representation_ptr->Surface()->
599               Transformed(local_location.Transformation()) ) ;
600           }
601           else {
602             surface_ptr = 
603               geometric_representation_ptr->Surface() ;
604           }
605           not_done = 0 ;
606       }
607       curve_on_surface_index += 1 ;
608     }
609     Geom2dAdaptor_Curve     AnAdaptor3dCurve2d (curve2d_ptr) ;
610     GeomAdaptor_Surface     AnAdaptor3dSurface (surface_ptr) ;
611     Handle(Geom2dAdaptor_HCurve) AnAdaptor3dCurve2dPtr =
612       new Geom2dAdaptor_HCurve(AnAdaptor3dCurve2d) ;
613     Handle(GeomAdaptor_HSurface) AnAdaptor3dSurfacePtr =
614       new GeomAdaptor_HSurface (AnAdaptor3dSurface) ;
615     curve_on_surface_reference.Load (AnAdaptor3dCurve2dPtr, AnAdaptor3dSurfacePtr);
616     a_sampler.Initialize(curve_on_surface_reference,
617       MinToleranceRequested * factor,
618       current_first,
619       current_last) ;
620   }
621   TColStd_Array1OfReal   sampling_parameters(1,a_sampler.NbPoints()) ;
622   for (ii = 1 ; ii <= a_sampler.NbPoints() ; ii++) {
623     sampling_parameters(ii) = a_sampler.Parameter(ii) ;
624   }
625   if (a_sampler.NbPoints() < min_sampling_points) {
626     GeomLib::DensifyArray1OfReal(min_sampling_points,
627       sampling_parameters,
628       parameters_ptr) ;
629   }
630   else if (a_sampler.NbPoints() > max_sampling_points) {
631     GeomLib::RemovePointsFromArray(max_sampling_points,
632       sampling_parameters,
633       parameters_ptr) ; 
634   }
635   else {
636     jj = 1 ;
637     parameters_ptr =
638       new TColStd_HArray1OfReal(1,sampling_parameters.Length()) ;
639     for (ii = sampling_parameters.Lower() ; ii <= sampling_parameters.Upper() ; ii++) {
640       parameters_ptr->ChangeArray1()(jj) =
641         sampling_parameters(ii) ;
642       jj +=1 ;
643     }
644   }
645
646   curve_index = 0 ;
647
648   while (second_iterator.More()) {
649     geometric_representation_ptr =
650       Handle(BRep_GCurve)::DownCast(second_iterator.Value());
651     if (! geometric_representation_ptr.IsNull() && 
652       curve_index != curve_on_surface_index) {
653         has_closed_curve =
654           has_curve = Standard_False ;
655         //      first = geometric_representation_ptr->First();
656         //      last =  geometric_representation_ptr->Last();
657         local_location = geometric_representation_ptr->Location() ;
658         if (geometric_representation_ptr->IsCurveOnSurface()) {
659           curve2d_ptr = geometric_representation_ptr->PCurve() ; 
660           has_curve = Standard_True ;
661         }
662         if (geometric_representation_ptr->IsCurveOnClosedSurface()) {
663           curve2d_ptr = geometric_representation_ptr->PCurve2() ;
664           has_closed_curve = Standard_True ;
665         }
666
667         if (has_curve ||
668           has_closed_curve) {
669             if (! local_location.IsIdentity() ) {
670               surface_ptr = Handle(Geom_Surface)::
671                 DownCast( geometric_representation_ptr->Surface()->
672                 Transformed(local_location.Transformation()) ) ;
673             }
674             else {
675               surface_ptr = 
676                 geometric_representation_ptr->Surface() ;
677             }
678             Geom2dAdaptor_Curve     an_adaptor_curve2d (curve2d_ptr) ;
679             GeomAdaptor_Surface     an_adaptor_surface(surface_ptr) ;
680             Handle(Geom2dAdaptor_HCurve) an_adaptor_curve2d_ptr =
681               new Geom2dAdaptor_HCurve(an_adaptor_curve2d) ;
682             Handle(GeomAdaptor_HSurface) an_adaptor_surface_ptr =
683               new GeomAdaptor_HSurface (an_adaptor_surface) ;
684             Adaptor3d_CurveOnSurface a_curve_on_surface(an_adaptor_curve2d_ptr,
685               an_adaptor_surface_ptr) ;
686
687             if (BRep_Tool::SameParameter(AnEdge)) {
688
689               GeomLib::EvalMaxParametricDistance(a_curve_on_surface,
690                 geom_reference_curve,
691                 MinToleranceRequested,
692                 parameters_ptr->Array1(),
693                 max_distance) ;
694             }
695             else if (geom_reference_curve_flag) {
696               GeomLib::EvalMaxDistanceAlongParameter(a_curve_on_surface,
697                 geom_reference_curve,
698                 MinToleranceRequested,
699                 parameters_ptr->Array1(),
700                 max_distance) ;
701             }
702             else {
703
704               GeomLib::EvalMaxDistanceAlongParameter(a_curve_on_surface,
705                 curve_on_surface_reference,
706                 MinToleranceRequested,
707                 parameters_ptr->Array1(),
708                 max_distance) ;
709             }
710             max_distance *= safe_factor ;
711             edge_tolerance = Max(max_distance, edge_tolerance) ;
712         }
713
714
715     }
716     curve_index += 1 ;
717     second_iterator.Next() ; 
718   }
719
720   TE->Tolerance(edge_tolerance);
721   return Standard_True ;
722
723 }
724 //=======================================================================
725 //function : UpdateEdgeTolerance
726 //purpose  : 
727 //=======================================================================
728
729 Standard_Boolean BRepLib::UpdateEdgeTolerance(const TopoDS_Shape& S,
730   const Standard_Real MinToleranceRequested,
731   const Standard_Real MaxToleranceToCheck) 
732 {
733   TopExp_Explorer ex(S,TopAbs_EDGE);
734   TopTools_MapOfShape  a_counter ;
735
736   Standard_Boolean     return_status = Standard_False,
737     local_flag ;
738
739   while (ex.More()) {
740     if (a_counter.Add(ex.Current())) {
741       local_flag =
742         BRepLib::UpdateEdgeTol(TopoDS::Edge(ex.Current()),
743         MinToleranceRequested,
744         MaxToleranceToCheck) ;
745       if (local_flag && ! return_status) {
746         return_status = Standard_True ;
747       }
748     }
749     ex.Next();
750   }
751   return return_status ;
752 }
753
754 //=======================================================================
755 //function : SetEdgeTol
756 //purpose  : 
757 //=======================================================================
758
759 static void SetEdgeTol(const TopoDS_Edge& E,
760   const TopoDS_Face& F)
761 {
762   TopLoc_Location L;
763   const Handle(Geom_Surface)& S = BRep_Tool::Surface(F,L);
764   TopLoc_Location l = L.Predivided(E.Location());
765
766   const Handle(BRep_TEdge)& TE = *((Handle(BRep_TEdge)*)&E.TShape());
767   BRep_ListIteratorOfListOfCurveRepresentation itcr(TE->ChangeCurves());
768
769   while (itcr.More()) {
770     const Handle(BRep_CurveRepresentation)& cr = itcr.Value();
771     if(cr->IsCurveOnSurface(S,l)) return;
772     itcr.Next();
773   }
774
775   Handle(Geom_Plane) GP;
776   Handle(Geom_RectangularTrimmedSurface) GRTS;
777   GRTS = Handle(Geom_RectangularTrimmedSurface)::DownCast(S);
778   if(!GRTS.IsNull())
779     GP = Handle(Geom_Plane)::DownCast(GRTS->BasisSurface());
780   else
781     GP = Handle(Geom_Plane)::DownCast(S);
782
783   Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve();
784   Handle(GeomAdaptor_HSurface) HS = new GeomAdaptor_HSurface();
785
786   TopLoc_Location LC;
787   Standard_Real First, Last;
788   GeomAdaptor_Curve& GAC = HC->ChangeCurve();
789   GAC.Load(BRep_Tool::Curve(E,LC,First,Last));
790   LC = L.Predivided(LC);
791
792   if (!LC.IsIdentity()) {
793     GP = Handle(Geom_Plane)::DownCast(
794       GP->Transformed(LC.Transformation()));
795   }
796   GeomAdaptor_Surface& GAS = HS->ChangeSurface();
797   GAS.Load(GP);
798
799   ProjLib_ProjectedCurve Proj(HS,HC);
800   Handle(Geom2d_Curve) pc = Geom2dAdaptor::MakeCurve(Proj);
801
802   gp_Pln pln = GAS.Plane();
803   Standard_Real d2 = 0.;
804   Standard_Integer nn = 23;
805   Standard_Real unsurnn = 1./nn;
806   for(Standard_Integer i = 0; i <= nn; i++){
807     Standard_Real t = unsurnn*i;
808     Standard_Real u = First*(1.-t) + Last*t;
809     gp_Pnt Pc3d = HC->Value(u);
810     gp_Pnt2d p2d = pc->Value(u);
811     gp_Pnt Pcons = ElSLib::Value(p2d.X(),p2d.Y(),pln);
812     Standard_Real eps = Max(Pc3d.XYZ().SquareModulus(), Pcons.XYZ().SquareModulus());
813     eps = Epsilon(eps);
814     Standard_Real temp = Pc3d.SquareDistance(Pcons);
815     if(temp <= eps)
816     {
817       temp = 0.;
818     }
819     if(temp > d2) d2 = temp;
820   }
821   d2 = 1.5*sqrt(d2);
822   TE->UpdateTolerance(d2);
823 }
824
825 //=======================================================================
826 //function : SameParameter
827 //purpose  : 
828 //=======================================================================
829 void  BRepLib::SameParameter(const TopoDS_Shape& S,
830   const Standard_Real Tolerance,
831   const Standard_Boolean forced) 
832 {
833   TopExp_Explorer ex(S,TopAbs_EDGE);
834   TopTools_MapOfShape  Done;
835   BRep_Builder brB;
836
837   while (ex.More()) {
838     if (Done.Add(ex.Current())) {
839       if (forced) {
840         brB.SameRange(TopoDS::Edge(ex.Current()), Standard_False);
841         brB.SameParameter(TopoDS::Edge(ex.Current()), Standard_False);
842       }
843       BRepLib::SameParameter(TopoDS::Edge(ex.Current()),Tolerance);
844     }
845     ex.Next();
846   }
847
848   Done.Clear();
849   BRepAdaptor_Surface BS;
850   for(ex.Init(S,TopAbs_FACE); ex.More(); ex.Next()){
851     const TopoDS_Face& curface = TopoDS::Face(ex.Current());
852     if(!Done.Add(curface)) continue;
853     BS.Initialize(curface);
854     if(BS.GetType() != GeomAbs_Plane) continue;
855     TopExp_Explorer ex2;
856     for(ex2.Init(curface,TopAbs_EDGE); ex2.More(); ex2.Next()){
857       const TopoDS_Edge& E = TopoDS::Edge(ex2.Current());
858       SetEdgeTol(E,curface);
859     }
860   }
861   BRepLib::UpdateTolerances(S);
862 }
863
864 //================================================================
865 //function : SameParameter
866 //WARNING  : New spec DUB LBO 9/9/97.
867 //  Recode in the edge the best tolerance found, 
868 //  for vertex extremities it is required to find something else
869 //================================================================
870 static Standard_Boolean EvalTol(const Handle(Geom2d_Curve)& pc,
871   const Handle(Geom_Surface)& s,
872   const GeomAdaptor_Curve&    gac,
873   const Standard_Real         tol,
874   Standard_Real&              tolbail)
875 {
876   Standard_Integer ok = 0;
877   Standard_Real f = gac.FirstParameter();
878   Standard_Real l = gac.LastParameter();
879   Extrema_LocateExtPC Projector;
880   Projector.Initialize(gac,f,l,tol);
881   Standard_Real u,v;
882   gp_Pnt p;
883   tolbail = tol;
884   for(Standard_Integer i = 1; i <= 5; i++){
885     Standard_Real t = i/6.;
886     t = (1.-t) * f + t * l;
887     pc->Value(t).Coord(u,v);
888     p = s->Value(u,v);
889     Projector.Perform(p,t);
890     if (Projector.IsDone()) {
891       Standard_Real dist2 = Projector.SquareDistance();
892       if(dist2 > tolbail * tolbail) tolbail = sqrt (dist2);
893       ok++;
894     }
895   }
896   return (ok > 2);
897 }
898
899 static Standard_Real ComputeTol(const Handle(Adaptor3d_HCurve)& c3d,
900   const Handle(Adaptor2d_HCurve2d)& c2d,
901   const Handle(Adaptor3d_HSurface)& surf,
902   const Standard_Integer        nbp)
903
904 {
905
906   TColStd_Array1OfReal dist(1,nbp+10);
907   dist.Init(-1.);
908
909   //Adaptor3d_CurveOnSurface  cons(c2d,surf);
910   Standard_Real uf = surf->FirstUParameter(), ul = surf->LastUParameter(),
911                 vf = surf->FirstVParameter(), vl = surf->LastVParameter();
912   Standard_Real du = 0.01 * (ul - uf), dv = 0.01 * (vl - vf);
913   Standard_Boolean isUPeriodic = surf->IsUPeriodic(), isVPeriodic = surf->IsVPeriodic();
914   Standard_Real DSdu = 1./surf->UResolution(1.), DSdv = 1./surf->VResolution(1.);
915   Standard_Real d2 = 0.;
916   Standard_Real first = c3d->FirstParameter();
917   Standard_Real last  = c3d->LastParameter();
918   Standard_Real dapp = -1.;
919   Standard_Integer i = 0;
920   for(i = 0; i <= nbp; i++){
921     const Standard_Real t = IntToReal(i)/IntToReal(nbp);
922     const Standard_Real u = first*(1.-t) + last*t;
923     gp_Pnt Pc3d = c3d->Value(u);
924     gp_Pnt2d Puv = c2d->Value(u);
925     if(!isUPeriodic)
926     {
927       if(Puv.X() < uf - du)
928       {
929         dapp = Max(dapp, DSdu * (uf - Puv.X()));
930         continue;
931       }
932       else if(Puv.X() > ul + du)
933       {
934         dapp = Max(dapp, DSdu * (Puv.X() - ul));
935         continue;
936       }
937     }
938     if(!isVPeriodic)
939     {
940       if(Puv.Y() < vf - dv)
941       {
942         dapp = Max(dapp, DSdv * (vf - Puv.Y()));
943         continue;
944       }
945       else if(Puv.Y() > vl + dv)
946       {
947         dapp = Max(dapp, DSdv * (Puv.Y() - vl));
948         continue;
949       }
950     }
951     gp_Pnt Pcons = surf->Value(Puv.X(), Puv.Y());
952     if (Precision::IsInfinite(Pcons.X()) ||
953         Precision::IsInfinite(Pcons.Y()) ||
954         Precision::IsInfinite(Pcons.Z()))
955     {
956         d2=Precision::Infinite();
957         break;
958     }
959     Standard_Real temp = Pc3d.SquareDistance(Pcons);
960
961
962     dist(i+1) = temp;
963
964
965     if(temp > d2) d2 = temp;
966   }
967
968   if(Precision::IsInfinite(d2))
969   {
970     return d2;
971   }
972
973   d2 = Sqrt(d2);
974   if(dapp > d2)
975   {
976     return dapp;
977   }
978
979   Standard_Boolean ana = Standard_False;
980   Standard_Real D2 = 0;
981   Standard_Integer N1 = 0;
982   Standard_Integer N2 = 0;
983   Standard_Integer N3 = 0;
984
985   for( i = 1; i<= nbp+10; i++)
986     if( dist(i) > 0 ) {
987       if( dist(i) < 1.0 ) N1++;
988       else N2++;
989     }
990
991     if( N1 > N2 && N2 != 0 ) N3 = 100*N2/(N1+N2);
992     if( N3 < 10 && N3 != 0 ) {
993       ana = Standard_True;
994       for( i = 1; i<= nbp+10; i++)
995         if( dist(i) > 0 && dist(i) < 1.0 )
996           if( dist(i) > D2 ) D2 = dist(i);
997     }
998
999     //d2 = 1.5*sqrt(d2);
1000     d2 = (!ana) ? 1.5 * d2 : 1.5*sqrt(D2);
1001     if(d2<1.e-7) d2 = 1.e-7;
1002
1003     return d2;
1004 }
1005
1006 void BRepLib::SameParameter(const TopoDS_Edge&  AnEdge,
1007   const Standard_Real Tolerance)
1008 {
1009   if (BRep_Tool::SameParameter(AnEdge)) return;
1010
1011   const Standard_Integer NCONTROL = 22;
1012
1013   Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve();
1014   Handle(Geom2dAdaptor_HCurve) HC2d = new Geom2dAdaptor_HCurve();
1015   Handle(GeomAdaptor_HSurface) HS = new GeomAdaptor_HSurface();
1016   GeomAdaptor_Curve& GAC = HC->ChangeCurve();
1017   Geom2dAdaptor_Curve& GAC2d = HC2d->ChangeCurve2d();
1018   GeomAdaptor_Surface& GAS = HS->ChangeSurface();
1019   Standard_Real f3d =0.,l3d =0.;
1020   TopLoc_Location L3d;
1021   Handle(Geom_Curve) C3d;
1022
1023   const Handle(BRep_TEdge)& TE = *((Handle(BRep_TEdge)*) &AnEdge.TShape());
1024   BRep_ListOfCurveRepresentation& CList = TE->ChangeCurves();
1025   BRep_ListIteratorOfListOfCurveRepresentation It(CList);
1026
1027   Standard_Boolean NotDone = Standard_True;
1028
1029   while (NotDone && It.More()) {
1030     Handle(BRep_GCurve) GCurve = Handle(BRep_GCurve)::DownCast(It.Value());
1031     if (!GCurve.IsNull() && GCurve->IsCurve3D()) {
1032       C3d = GCurve->Curve3D() ;
1033       f3d = GCurve->First();
1034       l3d = GCurve->Last();
1035       L3d = GCurve->Location() ;
1036       NotDone = Standard_False;
1037     } 
1038     It.Next() ;
1039   }
1040
1041   if(C3d.IsNull()) return;
1042
1043   // modified by NIZHNY-OCC486  Tue Aug 27 17:15:13 2002 :
1044   Standard_Boolean m_TrimmedPeriodical = Standard_False;
1045   Handle(Standard_Type) TheType = C3d->DynamicType();
1046   if( TheType == STANDARD_TYPE(Geom_TrimmedCurve))
1047   {
1048     Handle(Geom_Curve) gtC (Handle(Geom_TrimmedCurve)::DownCast (C3d)->BasisCurve());
1049     m_TrimmedPeriodical = gtC->IsPeriodic();
1050   }
1051   // modified by NIZHNY-OCC486  Tue Aug 27 17:15:17 2002 .
1052
1053   BRep_Builder B;
1054   if(!C3d->IsPeriodic()) {
1055     Standard_Real Udeb = C3d->FirstParameter();
1056     Standard_Real Ufin = C3d->LastParameter();
1057     // modified by NIZHNY-OCC486  Tue Aug 27 17:17:14 2002 :
1058     //if (Udeb > f3d) f3d = Udeb;
1059     //if (l3d > Ufin) l3d = Ufin;
1060     if(!m_TrimmedPeriodical)
1061     {
1062       if (Udeb > f3d) f3d = Udeb;
1063       if (l3d > Ufin) l3d = Ufin;
1064     }
1065     // modified by NIZHNY-OCC486  Tue Aug 27 17:17:55 2002 .
1066   }
1067   if(!L3d.IsIdentity()){
1068     C3d = Handle(Geom_Curve)::DownCast(C3d->Transformed(L3d.Transformation()));
1069   }
1070   GAC.Load(C3d,f3d,l3d);
1071
1072   Standard_Boolean IsSameP = 1;
1073   Standard_Real maxdist = 0.;
1074
1075   //  Modified by skv - Thu Jun  3 12:39:19 2004 OCC5898 Begin
1076   Standard_Real anEdgeTol = BRep_Tool::Tolerance(AnEdge);
1077   //  Modified by skv - Thu Jun  3 12:39:20 2004 OCC5898 End
1078   Standard_Boolean SameRange = BRep_Tool::SameRange(AnEdge);
1079   Standard_Boolean YaPCu = Standard_False;
1080   const Standard_Real BigError = 1.e10;
1081   It.Initialize(CList);
1082
1083   while (It.More()) {
1084     Standard_Boolean isANA = Standard_False;
1085     Standard_Boolean isBSP = Standard_False;
1086     Handle(BRep_GCurve) GCurve = Handle(BRep_GCurve)::DownCast(It.Value());
1087     Handle(Geom2d_Curve) PC[2];
1088     Handle(Geom_Surface) S;
1089     if (!GCurve.IsNull() && GCurve->IsCurveOnSurface()) {
1090       YaPCu = Standard_True;
1091       PC[0] = GCurve->PCurve();
1092       TopLoc_Location PCLoc = GCurve->Location();
1093       S = GCurve->Surface();
1094       if (!PCLoc.IsIdentity() ) {
1095         S = Handle(Geom_Surface)::DownCast(S->Transformed(PCLoc.Transformation()));
1096       }
1097
1098       GAS.Load(S);
1099       if (GCurve->IsCurveOnClosedSurface()) {
1100         PC[1] = GCurve->PCurve2();
1101       }
1102
1103       // Eval tol2d to compute SameRange
1104       Standard_Real UResol = Max(GAS.UResolution(Tolerance), Precision::PConfusion());
1105       Standard_Real VResol = Max(GAS.VResolution(Tolerance), Precision::PConfusion());
1106       Standard_Real Tol2d  = Min(UResol, VResol);
1107       for(Standard_Integer i = 0; i < 2; i++){
1108         Handle(Geom2d_Curve) curPC = PC[i];
1109         Standard_Boolean updatepc = 0;
1110         if(curPC.IsNull()) break;
1111         if(!SameRange){
1112           GeomLib::SameRange(Tol2d,
1113             PC[i],GCurve->First(),GCurve->Last(),
1114             f3d,l3d,curPC);
1115
1116           updatepc = (curPC != PC[i]);
1117
1118         }
1119         Standard_Boolean goodpc = 1;
1120         GAC2d.Load(curPC,f3d,l3d);
1121
1122         Standard_Real error = ComputeTol(HC, HC2d, HS, NCONTROL);
1123
1124         if(error > BigError)
1125         {
1126           maxdist = error;
1127           break;
1128         }
1129
1130         if(GAC2d.GetType() == GeomAbs_BSplineCurve && 
1131           GAC2d.Continuity() == GeomAbs_C0) {
1132             Handle(Geom2d_BSplineCurve) bs2d = GAC2d.BSpline();
1133             Handle(Geom2d_BSplineCurve) bs2dsov = bs2d;
1134             Standard_Real fC0 = bs2d->FirstParameter(), lC0 = bs2d->LastParameter();
1135             Standard_Boolean repar = Standard_True;
1136             gp_Pnt2d OriginPoint;
1137             bs2d->D0(fC0, OriginPoint);
1138             Geom2dConvert::C0BSplineToC1BSplineCurve(bs2d, Tol2d);
1139             isBSP = Standard_True; 
1140
1141             if(bs2d->IsPeriodic()) { // -------- IFV, Jan 2000
1142               gp_Pnt2d NewOriginPoint;
1143               bs2d->D0(bs2d->FirstParameter(), NewOriginPoint);
1144               if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1145                 Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) {
1146
1147                   TColStd_Array1OfReal Knotbs2d (1, bs2d->NbKnots());
1148                   bs2d->Knots(Knotbs2d);
1149
1150                   for(Standard_Integer Index = 1; Index <= bs2d->NbKnots(); Index++) {
1151                     bs2d->D0(Knotbs2d(Index), NewOriginPoint);
1152                     if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1153                       Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) continue;
1154
1155                     bs2d->SetOrigin(Index);
1156                     break;
1157                   }
1158               }
1159             }
1160
1161             if(bs2d->Continuity() == GeomAbs_C0) {
1162               Standard_Real tolbail;
1163               if(EvalTol(curPC,S,GAC,Tolerance,tolbail)){
1164                 bs2d = bs2dsov;
1165                 Standard_Real UResbail = GAS.UResolution(tolbail);
1166                 Standard_Real VResbail = GAS.VResolution(tolbail);
1167                 Standard_Real Tol2dbail  = Min(UResbail,VResbail);
1168                 bs2d->D0(bs2d->FirstParameter(), OriginPoint); 
1169
1170                 Standard_Integer nbp = bs2d->NbPoles();
1171                 TColgp_Array1OfPnt2d poles(1,nbp);
1172                 bs2d->Poles(poles);
1173                 gp_Pnt2d p = poles(1), p1;
1174                 Standard_Real d = Precision::Infinite();
1175                 for(Standard_Integer ip = 2; ip <= nbp; ip++) {
1176                   p1 = poles(ip);
1177                   d = Min(d,p.SquareDistance(p1));
1178                   p = p1;
1179                 }
1180                 d = sqrt(d)*.1;
1181
1182                 Tol2dbail = Max(Min(Tol2dbail,d),Tol2d);
1183
1184                 Geom2dConvert::C0BSplineToC1BSplineCurve(bs2d,Tol2dbail);
1185
1186                 if(bs2d->IsPeriodic()) { // -------- IFV, Jan 2000
1187                   gp_Pnt2d NewOriginPoint;
1188                   bs2d->D0(bs2d->FirstParameter(), NewOriginPoint);
1189                   if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1190                     Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) {
1191
1192                       TColStd_Array1OfReal Knotbs2d (1, bs2d->NbKnots());
1193                       bs2d->Knots(Knotbs2d);
1194
1195                       for(Standard_Integer Index = 1; Index <= bs2d->NbKnots(); Index++) {
1196                         bs2d->D0(Knotbs2d(Index), NewOriginPoint);
1197                         if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1198                           Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) continue;
1199
1200                         bs2d->SetOrigin(Index);
1201                         break;
1202                       }
1203                   }
1204                 }
1205
1206
1207                 if(bs2d->Continuity() == GeomAbs_C0) {
1208                   goodpc = 1;
1209                   bs2d = bs2dsov;
1210                   repar = Standard_False;
1211                 }
1212               }
1213               else goodpc = 0;
1214             }
1215
1216             if(goodpc){
1217               if(repar) {
1218                 Standard_Integer NbKnots = bs2d->NbKnots();
1219                 TColStd_Array1OfReal Knots(1,NbKnots);
1220                 bs2d->Knots(Knots);
1221                 //          BSplCLib::Reparametrize(f3d,l3d,Knots);
1222                 BSplCLib::Reparametrize(fC0,lC0,Knots);
1223                 bs2d->SetKnots(Knots);
1224                 GAC2d.Load(bs2d,f3d,l3d);
1225                 curPC = bs2d;
1226                 Standard_Boolean updatepcsov = updatepc;
1227                 updatepc = Standard_True;
1228
1229                 Standard_Real error1 = ComputeTol(HC, HC2d, HS, NCONTROL);
1230                 if(error1 > error) {
1231                   bs2d = bs2dsov;
1232                   GAC2d.Load(bs2d,f3d,l3d);
1233                   curPC = bs2d;
1234                   updatepc = updatepcsov;
1235                   isANA = Standard_True;
1236                 }
1237                 else {
1238                   error = error1;
1239                 }
1240               }
1241
1242               //check, if new BSpline "good" or not --------- IFV, Jan of 2000
1243               GeomAbs_Shape cont = bs2d->Continuity();
1244               Standard_Boolean IsBad = Standard_False;
1245
1246               if(cont > GeomAbs_C0 && error > Max(1.e-3,Tolerance)) {
1247                 Standard_Integer NbKnots = bs2d->NbKnots();
1248                 TColStd_Array1OfReal Knots(1,NbKnots);
1249                 bs2d->Knots(Knots);
1250                 Standard_Real critratio = 10.; 
1251                 Standard_Real dtprev = Knots(2) - Knots(1), dtratio = 1.;
1252                 Standard_Real dtmin = dtprev;
1253                 Standard_Real dtcur;
1254                 for(Standard_Integer j = 2; j < NbKnots; j++) {
1255                   dtcur = Knots(j+1) - Knots(j);
1256                   dtmin = Min(dtmin, dtcur);
1257
1258                   if(IsBad) continue;
1259
1260                   if(dtcur > dtprev) dtratio = dtcur/dtprev;
1261                   else dtratio = dtprev/dtcur;
1262                   if(dtratio > critratio) {IsBad = Standard_True;}
1263                   dtprev = dtcur;
1264
1265                 }
1266                 if(IsBad) {
1267                   // To avoid failures in Approx_CurvilinearParameter 
1268                   bs2d->Resolution(Max(1.e-3,Tolerance), dtcur);
1269                   if(dtmin < dtcur) IsBad = Standard_False;
1270                 }
1271               }
1272
1273
1274               if(IsBad ) { //if BSpline "bad", try to reparametrize it
1275                 // by its curve length
1276
1277                 //            GeomAbs_Shape cont = bs2d->Continuity();
1278                 if(cont > GeomAbs_C2) cont = GeomAbs_C2;
1279                 Standard_Integer maxdeg = bs2d->Degree();
1280                 if(maxdeg == 1) maxdeg = 14;
1281                 Approx_CurvilinearParameter AppCurPar(HC2d, HS, Max(1.e-3,Tolerance),
1282                   cont, maxdeg, 10);
1283                 if(AppCurPar.IsDone() || AppCurPar.HasResult()) {
1284                   bs2d = AppCurPar.Curve2d1();
1285                   GAC2d.Load(bs2d,f3d,l3d);
1286                   curPC = bs2d;
1287
1288                   if(Abs(bs2d->FirstParameter() - fC0) > Tol2d ||
1289                     Abs(bs2d->LastParameter() - lC0) > Tol2d    ) {
1290                       Standard_Integer NbKnots = bs2d->NbKnots();
1291                       TColStd_Array1OfReal Knots(1,NbKnots);
1292                       bs2d->Knots(Knots);
1293                       //                  BSplCLib::Reparametrize(f3d,l3d,Knots);
1294                       BSplCLib::Reparametrize(fC0,lC0,Knots);
1295                       bs2d->SetKnots(Knots);
1296                       GAC2d.Load(bs2d,f3d,l3d);
1297                       curPC = bs2d;
1298
1299                   }
1300                 }
1301               }
1302
1303
1304             }
1305         }
1306
1307
1308         if(goodpc){
1309           //      Approx_SameParameter SameP(HC,HC2d,HS,Tolerance);
1310           Standard_Real aTol = (isANA && isBSP) ? 1.e-7 : Tolerance;
1311           const Handle(Adaptor3d_HCurve)& aHCurv = HC; // to avoid ambiguity
1312           const Handle(Adaptor2d_HCurve2d)& aHCurv2d = HC2d; // to avoid ambiguity
1313           Approx_SameParameter SameP(aHCurv,aHCurv2d,HS,aTol);
1314
1315           if (SameP.IsSameParameter()) {
1316             maxdist = Max(maxdist,SameP.TolReached());
1317             if(updatepc){
1318               if (i == 0) GCurve->PCurve(curPC);
1319               else GCurve->PCurve2(curPC);
1320             }
1321           }
1322           else if (SameP.IsDone()) {
1323             Standard_Real tolreached = SameP.TolReached();
1324             if(tolreached <= error) {
1325               curPC = SameP.Curve2d();
1326               updatepc = Standard_True;
1327               maxdist = Max(maxdist,tolreached);
1328             }
1329             else {
1330               maxdist = Max(maxdist, error);
1331             }
1332             if(updatepc){
1333               if (i == 0) GCurve->PCurve(curPC);
1334               else GCurve->PCurve2(curPC);
1335             }
1336           }
1337           else
1338           {
1339             //Approx_SameParameter has failed.
1340             //Consequently, the situation might be,
1341             //when 3D and 2D-curve do not have same-range.
1342             GeomLib::SameRange( Tol2d, PC[i], 
1343                                 GCurve->First(), GCurve->Last(),
1344                                 f3d,l3d,curPC);
1345
1346             if (i == 0) GCurve->PCurve(curPC);
1347             else GCurve->PCurve2(curPC);
1348
1349             IsSameP = 0;
1350           }
1351
1352         }
1353         else IsSameP = 0;
1354
1355         //  Modified by skv - Thu Jun  3 12:39:19 2004 OCC5898 Begin
1356         if (!IsSameP) {
1357           if (anEdgeTol >= error) {
1358             maxdist = Max(maxdist, anEdgeTol);
1359             IsSameP = Standard_True;
1360           }
1361         }
1362         //  Modified by skv - Thu Jun  3 12:39:20 2004 OCC5898 End
1363       }
1364     }
1365     It.Next() ;
1366   }
1367   B.Range(AnEdge,f3d,l3d);
1368   B.SameRange(AnEdge,Standard_True);
1369   if ( IsSameP) {
1370     // Reduce eventually the tolerance of the edge, as
1371     // all its representations are processed (except for some associated
1372     // to planes and not stored in the edge !) 
1373     // The same cannot be done with vertices that cannot be enlarged 
1374     // or left as is.
1375     if (YaPCu) {
1376       // Avoid setting too small tolerances.
1377       maxdist = Max(maxdist,Precision::Confusion());
1378       TopoDS_Vertex V1,V2;
1379       TopExp::Vertices(AnEdge,V1,V2);
1380       if (!V1.IsNull())
1381         B.UpdateVertex(V1,maxdist);
1382       if (!V2.IsNull())
1383         B.UpdateVertex(V2,maxdist);
1384       TE->Modified(Standard_True);
1385       TE->Tolerance(maxdist);
1386     }
1387     B.SameParameter(AnEdge,Standard_True);
1388   }
1389 }
1390
1391 //=======================================================================
1392 //function : UpdateTolerances
1393 //purpose  : 
1394 //=======================================================================
1395 void  BRepLib::UpdateTolerances(const TopoDS_Shape& aShape,
1396   const Standard_Boolean verifyTolerance) 
1397 {
1398
1399   // Harmonize tolerances
1400   // with rule Tolerance(VERTEX)>=Tolerance(EDGE)>=Tolerance(FACE)
1401   BRep_Builder B;
1402   Standard_Real tol=0;
1403   if (verifyTolerance) {
1404     // Set tolerance to its minimum value
1405     Handle(Geom_Surface) S;
1406     TopLoc_Location l;
1407     TopExp_Explorer ex;
1408     Bnd_Box aB;
1409     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, dMax;
1410     for (ex.Init(aShape, TopAbs_FACE); ex.More(); ex.Next()) {
1411       const TopoDS_Face& curf=TopoDS::Face(ex.Current());
1412       S = BRep_Tool::Surface(curf, l);
1413       if (!S.IsNull()) {
1414         aB.SetVoid();
1415         BRepBndLib::Add(curf,aB);
1416         if (S->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
1417           S = Handle(Geom_RectangularTrimmedSurface)::DownCast (S)->BasisSurface();
1418         }
1419         GeomAdaptor_Surface AS(S);
1420         switch (AS.GetType()) {
1421         case GeomAbs_Plane: 
1422         case GeomAbs_Cylinder: 
1423         case GeomAbs_Cone: 
1424           {
1425             tol=Precision::Confusion();
1426             break;
1427           }
1428         case GeomAbs_Sphere: 
1429         case GeomAbs_Torus: 
1430           {
1431             tol=Precision::Confusion()*2;
1432             break;
1433           }
1434         default:
1435           tol=Precision::Confusion()*4;
1436         }
1437         if (!aB.IsWhole()) {
1438           aB.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
1439           dMax=1.;
1440           if (!aB.IsOpenXmin() && !aB.IsOpenXmax()) dMax=aXmax-aXmin;
1441           if (!aB.IsOpenYmin() && !aB.IsOpenYmax()) aYmin=aYmax-aYmin;
1442           if (!aB.IsOpenZmin() && !aB.IsOpenZmax()) aZmin=aZmax-aZmin;
1443           if (aYmin>dMax) dMax=aYmin;
1444           if (aZmin>dMax) dMax=aZmin;
1445           tol=tol*dMax;
1446           // Do not process tolerances > 1.
1447           if (tol>1.) tol=0.99;
1448         }
1449         const Handle(BRep_TFace)& Tf = *((Handle(BRep_TFace)*)&curf.TShape());
1450         Tf->Tolerance(tol);
1451       }
1452     }
1453   }
1454
1455   //Process edges
1456   TopTools_IndexedDataMapOfShapeListOfShape parents;
1457   TopExp::MapShapesAndAncestors(aShape, TopAbs_EDGE, TopAbs_FACE, parents);
1458   TopTools_ListIteratorOfListOfShape lConx;
1459   Standard_Integer iCur;
1460   for (iCur=1; iCur<=parents.Extent(); iCur++) {
1461     tol=0;
1462     for (lConx.Initialize(parents(iCur)); lConx.More(); lConx.Next()) {
1463       tol=Max(tol, BRep_Tool::Tolerance(TopoDS::Face(lConx.Value())));
1464     }
1465     // Update can only increase tolerance, so if the edge has a greater
1466     //  tolerance than its faces it is not concerned
1467     B.UpdateEdge(TopoDS::Edge(parents.FindKey(iCur)), tol);
1468   }
1469
1470   //Vertices are processed
1471   const Standard_Real BigTol = 1.e10;
1472   parents.Clear();
1473   TopExp::MapShapesAndAncestors(aShape, TopAbs_VERTEX, TopAbs_EDGE, parents);
1474   TColStd_MapOfTransient Initialized;
1475   TopTools_MapOfShape Done;
1476   Standard_Integer nbV = parents.Extent();
1477   for (iCur=1; iCur<=nbV; iCur++) {
1478     tol=0;
1479     Done.Clear();
1480     const TopoDS_Vertex& V = TopoDS::Vertex(parents.FindKey(iCur));
1481     Bnd_Box box;
1482     box.Add(BRep_Tool::Pnt(V));
1483     gp_Pnt p3d;
1484     for (lConx.Initialize(parents(iCur)); lConx.More(); lConx.Next()) {
1485       const TopoDS_Edge& E = TopoDS::Edge(lConx.Value());
1486       if(!Done.Add(E)) continue;
1487       tol=Max(tol, BRep_Tool::Tolerance(E));
1488       if(tol > BigTol) continue;
1489       if(!BRep_Tool::SameRange(E)) continue;
1490       Standard_Real par = BRep_Tool::Parameter(V,E);
1491       Handle(BRep_TEdge)& TE = *((Handle(BRep_TEdge)*)&E.TShape());
1492       BRep_ListIteratorOfListOfCurveRepresentation itcr(TE->Curves());
1493       const TopLoc_Location& Eloc = E.Location();
1494       while (itcr.More()) {
1495         // For each CurveRepresentation, check the provided parameter
1496         const Handle(BRep_CurveRepresentation)& cr = itcr.Value();
1497         const TopLoc_Location& loc = cr->Location();
1498         TopLoc_Location L = (Eloc * loc);
1499         if (cr->IsCurve3D()) {
1500           const Handle(Geom_Curve)& C = cr->Curve3D();
1501           if (!C.IsNull()) { // edge non degenerated
1502             p3d = C->Value(par);
1503             p3d.Transform(L.Transformation());
1504             box.Add(p3d);
1505           }
1506         }
1507         else if (cr->IsCurveOnSurface()) {
1508           const Handle(Geom_Surface)& Su = cr->Surface();
1509           const Handle(Geom2d_Curve)& PC = cr->PCurve();
1510           Handle(Geom2d_Curve) PC2;
1511           if (cr->IsCurveOnClosedSurface()) {
1512             PC2 = cr->PCurve2();
1513           }
1514           gp_Pnt2d p2d = PC->Value(par);
1515           p3d = Su->Value(p2d.X(),p2d.Y());
1516           p3d.Transform(L.Transformation());
1517           box.Add(p3d);
1518           if (!PC2.IsNull()) {
1519             p2d = PC2->Value(par);
1520             p3d = Su->Value(p2d.X(),p2d.Y());
1521             p3d.Transform(L.Transformation());
1522             box.Add(p3d);
1523           }
1524         }
1525         itcr.Next();
1526       }
1527     }
1528     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
1529     box.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
1530     aXmax -= aXmin; aYmax -= aYmin; aZmax -= aZmin;
1531     tol = Max(tol,sqrt(aXmax*aXmax+aYmax*aYmax+aZmax*aZmax));
1532     tol += 2.*Epsilon(tol);
1533     if (verifyTolerance) {
1534       // ASet minimum value of the tolerance 
1535       // Attention to sharing of the vertex by other shapes
1536       const Handle(BRep_TVertex)& TV = *((Handle(BRep_TVertex)*)&V.TShape());
1537       if (Initialized.Add(TV)) 
1538         TV->Tolerance(tol);
1539       else 
1540         B.UpdateVertex(V, tol);
1541     }
1542     else {
1543       // Update can only increase tolerance, so if the edge has a greater
1544       //  tolerance than its faces it is not concerned
1545       B.UpdateVertex(V, tol);
1546     }
1547   }
1548 }
1549
1550 //=======================================================================
1551 //function : UpdateInnerTolerances
1552 //purpose  : 
1553 //=======================================================================
1554 void  BRepLib::UpdateInnerTolerances(const TopoDS_Shape& aShape)
1555 {
1556   TopTools_IndexedDataMapOfShapeListOfShape EFmap;
1557   TopExp::MapShapesAndAncestors(aShape, TopAbs_EDGE, TopAbs_FACE, EFmap);
1558   BRep_Builder BB;
1559   for (Standard_Integer i = 1; i <= EFmap.Extent(); i++)
1560   {
1561     TopoDS_Edge anEdge = TopoDS::Edge(EFmap.FindKey(i));
1562     TopoDS_Vertex V1, V2;
1563     TopExp::Vertices(anEdge, V1, V2);
1564     Standard_Real fpar, lpar;
1565     BRep_Tool::Range(anEdge, fpar, lpar);
1566     Standard_Real TolEdge = BRep_Tool::Tolerance(anEdge);
1567     gp_Pnt Pnt1, Pnt2;
1568     Handle(BRepAdaptor_HCurve) anHCurve = new BRepAdaptor_HCurve();
1569     anHCurve->ChangeCurve().Initialize(anEdge);
1570     if (!V1.IsNull())
1571       Pnt1 = BRep_Tool::Pnt(V1);
1572     if (!V2.IsNull())
1573       Pnt2 = BRep_Tool::Pnt(V2);
1574     
1575     if (!BRep_Tool::Degenerated(anEdge) &&
1576         EFmap(i).Extent() > 0)
1577     {
1578       NCollection_Sequence<Handle(Adaptor3d_HCurve)> theRep;
1579       theRep.Append(anHCurve);
1580       TopTools_ListIteratorOfListOfShape itl(EFmap(i));
1581       for (; itl.More(); itl.Next())
1582       {
1583         const TopoDS_Face& aFace = TopoDS::Face(itl.Value());
1584         Handle(BRepAdaptor_HCurve) anHCurvOnSurf = new BRepAdaptor_HCurve();
1585         anHCurvOnSurf->ChangeCurve().Initialize(anEdge, aFace);
1586         theRep.Append(anHCurvOnSurf);
1587       }
1588       
1589       const Standard_Integer NbSamples = (BRep_Tool::SameParameter(anEdge))? 23 : 2;
1590       Standard_Real delta = (lpar - fpar)/(NbSamples-1);
1591       Standard_Real MaxDist = 0.;
1592       for (Standard_Integer j = 2; j <= theRep.Length(); j++)
1593       {
1594         for (Standard_Integer k = 0; k <= NbSamples; k++)
1595         {
1596           Standard_Real ParamOnCenter = (k == NbSamples)? lpar :
1597             fpar + k*delta;
1598           gp_Pnt Center = theRep(1)->Value(ParamOnCenter);
1599           Standard_Real ParamOnCurve = (BRep_Tool::SameParameter(anEdge))? ParamOnCenter
1600             : ((k == 0)? theRep(j)->FirstParameter() : theRep(j)->LastParameter());
1601           gp_Pnt aPoint = theRep(j)->Value(ParamOnCurve);
1602           Standard_Real aDist = Center.Distance(aPoint);
1603           //aDist *= 1.1;
1604           aDist += 2.*Epsilon(aDist);
1605           if (aDist > MaxDist)
1606             MaxDist = aDist;
1607
1608           //Update tolerances of vertices
1609           if (k == 0 && !V1.IsNull())
1610           {
1611             Standard_Real aDist1 = Pnt1.Distance(aPoint);
1612             aDist1 += 2.*Epsilon(aDist1);
1613             BB.UpdateVertex(V1, aDist1);
1614           }
1615           if (k == NbSamples && !V2.IsNull())
1616           {
1617             Standard_Real aDist2 = Pnt2.Distance(aPoint);
1618             aDist2 += 2.*Epsilon(aDist2);
1619             BB.UpdateVertex(V2, aDist2);
1620           }
1621         }
1622       }
1623       BB.UpdateEdge(anEdge, MaxDist);
1624     }
1625     TolEdge = BRep_Tool::Tolerance(anEdge);
1626     if (!V1.IsNull())
1627     {
1628       gp_Pnt End1 = anHCurve->Value(fpar);
1629       Standard_Real dist1 = Pnt1.Distance(End1);
1630       dist1 += 2.*Epsilon(dist1);
1631       BB.UpdateVertex(V1, Max(dist1, TolEdge));
1632     }
1633     if (!V2.IsNull())
1634     {
1635       gp_Pnt End2 = anHCurve->Value(lpar);
1636       Standard_Real dist2 = Pnt2.Distance(End2);
1637       dist2 += 2.*Epsilon(dist2);
1638       BB.UpdateVertex(V2, Max(dist2, TolEdge));
1639     }
1640   }
1641 }
1642
1643 //=======================================================================
1644 //function : OrientClosedSolid
1645 //purpose  : 
1646 //=======================================================================
1647 Standard_Boolean BRepLib::OrientClosedSolid(TopoDS_Solid& solid) 
1648 {
1649   // Set material inside the solid
1650   BRepClass3d_SolidClassifier where(solid);
1651   where.PerformInfinitePoint(Precision::Confusion());
1652   if (where.State()==TopAbs_IN) {
1653     solid.Reverse();
1654   }
1655   else if (where.State()==TopAbs_ON || where.State()==TopAbs_UNKNOWN) 
1656     return Standard_False;
1657
1658   return Standard_True;
1659 }
1660
1661 // Structure for calculation of properties, necessary for decision about continuity
1662 class SurfaceProperties
1663 {
1664 public:
1665   SurfaceProperties(const Handle(Geom_Surface)& theSurface,
1666                     const gp_Trsf&              theSurfaceTrsf,
1667                     const Handle(Geom2d_Curve)& theCurve2D,
1668                     const Standard_Boolean      theReversed)
1669     : mySurfaceProps(theSurface, 2, Precision::Confusion()),
1670       mySurfaceTrsf(theSurfaceTrsf),
1671       myCurve2d(theCurve2D),
1672       myIsReversed(theReversed)
1673   {}
1674
1675   // Calculate derivatives on surface related to the point on curve
1676   void Calculate(const Standard_Real theParamOnCurve)
1677   {
1678     gp_Pnt2d aUV;
1679     myCurve2d->D1(theParamOnCurve, aUV, myCurveTangent);
1680     mySurfaceProps.SetParameters(aUV.X(), aUV.Y());
1681   }
1682
1683   // Returns point just calculated
1684   gp_Pnt Value() 
1685   { return mySurfaceProps.Value().Transformed(mySurfaceTrsf); }
1686
1687   // Calculate a derivative orthogonal to curve's tangent vector
1688   gp_Vec Derivative()
1689   {
1690     gp_Vec aDeriv;
1691     // direction orthogonal to tangent vector of the curve
1692     gp_Vec2d anOrtho(-myCurveTangent.Y(), myCurveTangent.X());
1693     Standard_Real aLen = anOrtho.Magnitude();
1694     if (aLen < Precision::Confusion())
1695       return aDeriv;
1696     anOrtho /= aLen;
1697     if (myIsReversed)
1698       anOrtho.Reverse();
1699
1700     aDeriv.SetLinearForm(anOrtho.X(), mySurfaceProps.D1U(),
1701                          anOrtho.Y(), mySurfaceProps.D1V());
1702     return aDeriv.Transformed(mySurfaceTrsf);
1703   }
1704
1705   // Calculate principal curvatures, which consist of minimal and maximal normal curvatures and
1706   // the directions on the tangent plane (principal direction) where the extremums are reached
1707   void Curvature(gp_Dir& thePrincipalDir1, Standard_Real& theCurvature1,
1708                  gp_Dir& thePrincipalDir2, Standard_Real& theCurvature2)
1709   {
1710     mySurfaceProps.CurvatureDirections(thePrincipalDir1, thePrincipalDir2);
1711     theCurvature1 = mySurfaceProps.MaxCurvature();
1712     theCurvature2 = mySurfaceProps.MinCurvature();
1713     if (myIsReversed)
1714     {
1715       theCurvature1 = -theCurvature1;
1716       theCurvature2 = -theCurvature2;
1717     }
1718     thePrincipalDir1.Transform(mySurfaceTrsf);
1719     thePrincipalDir2.Transform(mySurfaceTrsf);
1720   }
1721
1722 private:
1723   GeomLProp_SLProps    mySurfaceProps; // properties calculator
1724   gp_Trsf              mySurfaceTrsf;
1725   Handle(Geom2d_Curve) myCurve2d;
1726   Standard_Boolean     myIsReversed; // the face based on the surface is reversed
1727
1728   // tangent vector to Pcurve in UV
1729   gp_Vec2d myCurveTangent;
1730 };
1731
1732 //=======================================================================
1733 //function : tgtfaces
1734 //purpose  : check the angle at the border between two squares.
1735 //           Two shares should have a shared front edge.
1736 //=======================================================================
1737 static GeomAbs_Shape tgtfaces(const TopoDS_Edge& Ed,
1738                               const TopoDS_Face& F1,
1739                               const TopoDS_Face& F2,
1740                               const Standard_Real theAngleTol)
1741 {
1742   Standard_Boolean isSeam = F1.IsEqual(F2);
1743
1744   TopoDS_Edge E = Ed;
1745
1746   // Check if pcurves exist on both faces of edge
1747   Standard_Real aFirst,aLast;
1748   E.Orientation(TopAbs_FORWARD);
1749   Handle(Geom2d_Curve) aCurve1 = BRep_Tool::CurveOnSurface(E, F1, aFirst, aLast);
1750   if(aCurve1.IsNull())
1751     return GeomAbs_C0;
1752
1753   if (isSeam)
1754     E.Orientation(TopAbs_REVERSED);
1755   Handle(Geom2d_Curve) aCurve2 = BRep_Tool::CurveOnSurface(E, F2, aFirst, aLast);
1756   if(aCurve2.IsNull())
1757     return GeomAbs_C0;
1758
1759   TopLoc_Location aLoc1, aLoc2;
1760   Handle(Geom_Surface) aSurface1 = BRep_Tool::Surface(F1, aLoc1);
1761   const gp_Trsf& aSurf1Trsf = aLoc1.Transformation();
1762   Handle(Geom_Surface) aSurface2 = BRep_Tool::Surface(F2, aLoc2);
1763   const gp_Trsf& aSurf2Trsf = aLoc2.Transformation();
1764
1765   if (aSurface1->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
1766     aSurface1 = Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface1)->BasisSurface();
1767   if (aSurface2->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
1768     aSurface2 = Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface2)->BasisSurface();
1769
1770   // seam edge on elementary surface is always CN
1771   Standard_Boolean isElementary =
1772     (aSurface1->IsKind(STANDARD_TYPE(Geom_ElementarySurface)) &&
1773      aSurface2->IsKind(STANDARD_TYPE(Geom_ElementarySurface)));
1774   if (isSeam && isElementary)
1775   {
1776     return GeomAbs_CN;
1777   }
1778
1779   SurfaceProperties aSP1(aSurface1, aSurf1Trsf, aCurve1, F1.Orientation() == TopAbs_REVERSED);
1780   SurfaceProperties aSP2(aSurface2, aSurf2Trsf, aCurve2, F2.Orientation() == TopAbs_REVERSED);
1781
1782   Standard_Real f, l, eps;
1783   BRep_Tool::Range(E,f,l);
1784   Extrema_LocateExtPC ext;
1785   Handle(BRepAdaptor_HCurve) aHC2;
1786
1787   eps = (l - f)/100.;
1788   f += eps; // to avoid calculations on  
1789   l -= eps; // points of pointed squares.
1790
1791   const Standard_Real anAngleTol2 = theAngleTol * theAngleTol;
1792
1793   gp_Vec aDer1, aDer2;
1794   gp_Vec aNorm1;
1795   Standard_Real aSqLen1, aSqLen2;
1796   gp_Dir aCrvDir1[2], aCrvDir2[2];
1797   Standard_Real aCrvLen1[2], aCrvLen2[2];
1798
1799   GeomAbs_Shape aCont = (isElementary ? GeomAbs_CN : GeomAbs_C2);
1800   GeomAbs_Shape aCurCont;
1801   Standard_Real u;
1802   for (Standard_Integer i = 0; i <= 20 && aCont > GeomAbs_C0; i++)
1803   {
1804     // First suppose that this is sameParameter
1805     u = f + (l-f)*i/20;
1806
1807     // Check conditions for G1 and C1 continuity:
1808     // * calculate a derivative in tangent plane of each surface
1809     //   orthogonal to curve's tangent vector
1810     // * continuity is C1 if the vectors are equal
1811     // * continuity is G1 if the vectors are just parallel
1812     aCurCont = GeomAbs_C0;
1813
1814     aSP1.Calculate(u);
1815     aSP2.Calculate(u);
1816
1817     aDer1 = aSP1.Derivative();
1818     aSqLen1 = aDer1.SquareMagnitude();
1819     aDer2 = aSP2.Derivative();
1820     aSqLen2 = aDer2.SquareMagnitude();
1821     Standard_Boolean isSmoothSuspect = (aDer1.CrossSquareMagnitude(aDer2) <= anAngleTol2 * aSqLen1 * aSqLen2);
1822     if (!isSmoothSuspect)
1823     {
1824       // Refine by projection
1825       if (aHC2.IsNull())
1826       {
1827         // adaptor for pcurve on the second surface
1828         aHC2 = new BRepAdaptor_HCurve(BRepAdaptor_Curve(E, F2));
1829         ext.Initialize(aHC2->Curve(), f, l, Precision::PConfusion());
1830       }
1831       ext.Perform(aSP1.Value(), u);
1832       if (ext.IsDone() && ext.IsMin())
1833       {
1834         const Extrema_POnCurv& poc = ext.Point();
1835         aSP2.Calculate(poc.Parameter());
1836         aDer2 = aSP2.Derivative();
1837         aSqLen2 = aDer2.SquareMagnitude();
1838       }
1839       isSmoothSuspect = (aDer1.CrossSquareMagnitude(aDer2) <= anAngleTol2 * aSqLen1 * aSqLen2);
1840     }
1841     if (isSmoothSuspect)
1842     {
1843       aCurCont = GeomAbs_G1;
1844       if (Abs(Sqrt(aSqLen1) - Sqrt(aSqLen2)) < Precision::Confusion() &&
1845           aDer1.Dot(aDer2) > Precision::SquareConfusion()) // <= check vectors are codirectional
1846         aCurCont = GeomAbs_C1;
1847     }
1848     else
1849       return GeomAbs_C0;
1850
1851     if (aCont < GeomAbs_G2)
1852       continue; // no need further processing, because maximal continuity is less than G2
1853
1854     // Check conditions for G2 and C2 continuity:
1855     // * calculate principal curvatures on each surface
1856     // * continuity is C2 if directions of principal curvatures are equal on differenct surfaces
1857     // * continuity is G2 if directions of principal curvatures are just parallel
1858     //   and values of curvatures are the same
1859     aSP1.Curvature(aCrvDir1[0], aCrvLen1[0], aCrvDir1[1], aCrvLen1[1]);
1860     aSP2.Curvature(aCrvDir2[0], aCrvLen2[0], aCrvDir2[1], aCrvLen2[1]);
1861     for (Standard_Integer aStep = 0; aStep <= 1; ++aStep)
1862     {
1863       if (aCrvDir1[0].XYZ().CrossSquareMagnitude(aCrvDir2[aStep].XYZ()) <= Precision::SquareConfusion() &&
1864           Abs(aCrvLen1[0] - aCrvLen2[aStep]) < Precision::Confusion() &&
1865           aCrvDir1[1].XYZ().CrossSquareMagnitude(aCrvDir2[1 - aStep].XYZ()) <= Precision::SquareConfusion() &&
1866           Abs(aCrvLen1[1] - aCrvLen2[1 - aStep]) < Precision::Confusion())
1867       {
1868         if (aCurCont == GeomAbs_C1 &&
1869             aCrvDir1[0].Dot(aCrvDir2[aStep]) > Precision::Confusion() &&
1870             aCrvDir1[1].Dot(aCrvDir2[1 - aStep]) > Precision::Confusion())
1871           aCurCont = GeomAbs_C2;
1872         else
1873           aCurCont = GeomAbs_G2;
1874         break;
1875       }
1876     }
1877
1878     if (aCurCont < aCont)
1879       aCont = aCurCont;
1880   }
1881
1882   // according to the list of supported elementary surfaces,
1883   // if the continuity is C2, than it is totally CN
1884   if (isElementary && aCont == GeomAbs_C2)
1885     aCont = GeomAbs_CN;
1886   return aCont;
1887 }
1888
1889 //=======================================================================
1890 // function : EncodeRegularity
1891 // purpose  : Code the regularities on all edges of the shape, boundary of 
1892 //            two faces that do not have it.
1893 //            Takes into account that compound may consists of same solid
1894 //            placed with different transformations
1895 //=======================================================================
1896 static void EncodeRegularity(const TopoDS_Shape&        theShape,
1897                              const Standard_Real        theTolAng,
1898                              TopTools_MapOfShape&       theMap,
1899                              const TopTools_MapOfShape& theEdgesToEncode = TopTools_MapOfShape())
1900 {
1901   TopoDS_Shape aShape = theShape;
1902   TopLoc_Location aNullLoc;
1903   aShape.Location(aNullLoc); // nullify location
1904   if (!theMap.Add(aShape))
1905     return; // do not need to process shape twice
1906
1907   if (aShape.ShapeType() == TopAbs_COMPOUND ||
1908       aShape.ShapeType() == TopAbs_COMPSOLID)
1909   {
1910     for (TopoDS_Iterator it(aShape); it.More(); it.Next())
1911       EncodeRegularity(it.Value(), theTolAng, theMap, theEdgesToEncode);
1912     return;
1913   }
1914
1915   try {
1916     OCC_CATCH_SIGNALS
1917
1918     TopTools_IndexedDataMapOfShapeListOfShape M;
1919     TopExp::MapShapesAndAncestors(aShape, TopAbs_EDGE, TopAbs_FACE, M);
1920     TopTools_ListIteratorOfListOfShape It;
1921     TopExp_Explorer Ex;
1922     TopoDS_Face F1,F2;
1923     Standard_Boolean found;
1924     for (Standard_Integer i = 1; i <= M.Extent(); i++){
1925       TopoDS_Edge E = TopoDS::Edge(M.FindKey(i));
1926       if (!theEdgesToEncode.IsEmpty())
1927       {
1928         // process only the edges from the list to update their regularity
1929         TopoDS_Shape aPureEdge = E.Located(aNullLoc);
1930         aPureEdge.Orientation(TopAbs_FORWARD);
1931         if (!theEdgesToEncode.Contains(aPureEdge))
1932           continue;
1933       }
1934
1935       found = Standard_False;                                     
1936       F1.Nullify();
1937       for (It.Initialize(M.FindFromIndex(i)); It.More() && !found; It.Next()){
1938         if (F1.IsNull()) { F1 = TopoDS::Face(It.Value()); }
1939         else {
1940           const TopoDS_Face& aTmpF2 = TopoDS::Face(It.Value());
1941           if (!F1.IsSame(aTmpF2)){
1942             found = Standard_True;
1943             F2 = aTmpF2;
1944           }
1945         }
1946       }
1947       if (!found && !F1.IsNull()){//is it a sewing edge?
1948         TopAbs_Orientation orE = E.Orientation();
1949         TopoDS_Edge curE;
1950         for (Ex.Init(F1, TopAbs_EDGE); Ex.More() && !found; Ex.Next()){
1951           curE = TopoDS::Edge(Ex.Current());
1952           if (E.IsSame(curE) && orE != curE.Orientation()) {
1953             found = Standard_True;
1954             F2 = F1;
1955           }
1956         }
1957       }
1958       if (found)
1959         BRepLib::EncodeRegularity(E, F1, F2, theTolAng);
1960     }
1961   }
1962   catch (Standard_Failure const& anException) {
1963 #ifdef OCCT_DEBUG
1964     cout << "Warning: Exception in BRepLib::EncodeRegularity(): ";
1965     anException.Print(cout);
1966     cout << endl;
1967 #endif
1968     (void)anException;
1969   }
1970 }
1971
1972 //=======================================================================
1973 // function : EncodeRegularity
1974 // purpose  : code the regularities on all edges of the shape, boundary of 
1975 //            two faces that do not have it.
1976 //=======================================================================
1977
1978 void BRepLib::EncodeRegularity(const TopoDS_Shape& S,
1979   const Standard_Real TolAng)
1980 {
1981   TopTools_MapOfShape aMap;
1982   ::EncodeRegularity(S, TolAng, aMap);
1983 }
1984
1985 //=======================================================================
1986 // function : EncodeRegularity
1987 // purpose  : code the regularities on all edges in the list that do not 
1988 //            have it, and which are boundary of two faces on the shape.
1989 //=======================================================================
1990
1991 void BRepLib::EncodeRegularity(const TopoDS_Shape& S,
1992   const TopTools_ListOfShape& LE,
1993   const Standard_Real TolAng)
1994 {
1995   // Collect edges without location and orientation
1996   TopTools_MapOfShape aPureEdges;
1997   TopLoc_Location aNullLoc;
1998   TopTools_ListIteratorOfListOfShape anEdgeIt(LE);
1999   for (; anEdgeIt.More(); anEdgeIt.Next())
2000   {
2001     TopoDS_Shape anEdge = anEdgeIt.Value();
2002     anEdge.Location(aNullLoc);
2003     anEdge.Orientation(TopAbs_FORWARD);
2004     aPureEdges.Add(anEdge);
2005   }
2006
2007   TopTools_MapOfShape aMap;
2008   ::EncodeRegularity(S, TolAng, aMap, aPureEdges);
2009 }
2010
2011 //=======================================================================
2012 // function : EncodeRegularity
2013 // purpose  : code the regularity between 2 faces connected by edge 
2014 //=======================================================================
2015
2016 void BRepLib::EncodeRegularity(TopoDS_Edge& E,
2017   const TopoDS_Face& F1,
2018   const TopoDS_Face& F2,
2019   const Standard_Real TolAng)
2020 {
2021   BRep_Builder B;
2022   if(BRep_Tool::Continuity(E,F1,F2)<=GeomAbs_C0){
2023     try {
2024       GeomAbs_Shape aCont = tgtfaces(E, F1, F2, TolAng);
2025       B.Continuity(E,F1,F2,aCont);
2026       
2027     }
2028     catch(Standard_Failure)
2029     {
2030 #ifdef OCCT_DEBUG
2031       cout << "Failure: Exception in BRepLib::EncodeRegularity" << endl;
2032 #endif
2033     }
2034   }
2035 }
2036
2037 //=======================================================================
2038 // function : EnsureNormalConsistency
2039 // purpose  : Corrects the normals in Poly_Triangulation of faces.
2040 //              Returns TRUE if any correction is done.
2041 //=======================================================================
2042 Standard_Boolean BRepLib::
2043             EnsureNormalConsistency(const TopoDS_Shape& theShape,
2044                                     const Standard_Real theAngTol,
2045                                     const Standard_Boolean theForceComputeNormals)
2046 {
2047   const Standard_Real aThresDot = cos(theAngTol);
2048
2049   Standard_Boolean aRetVal = Standard_False, isNormalsFound = Standard_False;
2050
2051   // compute normals if they are absent
2052   TopExp_Explorer anExpFace(theShape,TopAbs_FACE);
2053   for (; anExpFace.More(); anExpFace.Next())
2054   {
2055     const TopoDS_Face& aFace = TopoDS::Face(anExpFace.Current());
2056     const Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
2057     if(aSurf.IsNull())
2058       continue;
2059     TopLoc_Location aLoc;
2060     const Handle(Poly_Triangulation)& aPT = BRep_Tool::Triangulation(aFace, aLoc);
2061     if(aPT.IsNull())
2062       continue;
2063     if (!theForceComputeNormals && aPT->HasNormals())
2064     {
2065       isNormalsFound = Standard_True;
2066       continue;
2067     }
2068
2069     GeomLProp_SLProps aSLP(aSurf, 2, Precision::Confusion());
2070     const Standard_Integer anArrDim = 3*aPT->NbNodes();
2071     Handle(TShort_HArray1OfShortReal) aNormArr = new TShort_HArray1OfShortReal(1, anArrDim);
2072     Standard_Integer anNormInd = aNormArr->Lower();
2073     for(Standard_Integer i = aPT->UVNodes().Lower(); i <= aPT->UVNodes().Upper(); i++)
2074     {
2075       const gp_Pnt2d &aP2d = aPT->UVNodes().Value(i);
2076       aSLP.SetParameters(aP2d.X(), aP2d.Y());
2077
2078       gp_XYZ aNorm(0.,0.,0.);
2079       if(!aSLP.IsNormalDefined())
2080       {
2081 #ifdef OCCT_DEBUG
2082         cout << "BRepLib::EnsureNormalConsistency(): Cannot find normal!" << endl;
2083 #endif
2084       }
2085       else
2086       {
2087         aNorm = aSLP.Normal().XYZ();
2088         if (aFace.Orientation() == TopAbs_REVERSED)
2089           aNorm.Reverse();
2090       }
2091       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.X());
2092       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.Y());
2093       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.Z());
2094     }
2095
2096     aRetVal = Standard_True;
2097     isNormalsFound = Standard_True;
2098     aPT->SetNormals(aNormArr);
2099   }
2100
2101   if(!isNormalsFound)
2102   {
2103     return aRetVal;
2104   }
2105
2106   // loop by edges
2107   TopTools_IndexedDataMapOfShapeListOfShape aMapEF;
2108   TopExp::MapShapesAndAncestors(theShape,TopAbs_EDGE,TopAbs_FACE,aMapEF);
2109   for(Standard_Integer anInd = 1; anInd <= aMapEF.Extent(); anInd++)
2110   {
2111     const TopoDS_Edge& anEdg = TopoDS::Edge(aMapEF.FindKey(anInd));
2112     const TopTools_ListOfShape& anEdgList = aMapEF.FindFromIndex(anInd);
2113     if (anEdgList.Extent() != 2)
2114       continue;
2115     TopTools_ListIteratorOfListOfShape anItF(anEdgList);
2116     const TopoDS_Face aFace1 = TopoDS::Face(anItF.Value());
2117     anItF.Next();
2118     const TopoDS_Face aFace2 = TopoDS::Face(anItF.Value());
2119     TopLoc_Location aLoc1, aLoc2;
2120     const Handle(Poly_Triangulation)& aPT1 = BRep_Tool::Triangulation(aFace1, aLoc1);
2121     const Handle(Poly_Triangulation)& aPT2 = BRep_Tool::Triangulation(aFace2, aLoc2);
2122
2123     if(aPT1.IsNull() || aPT2.IsNull())
2124       continue;
2125
2126     if(!aPT1->HasNormals() || !aPT2->HasNormals())
2127       continue;
2128
2129     const Handle(Poly_PolygonOnTriangulation)& aPTEF1 = 
2130                                 BRep_Tool::PolygonOnTriangulation(anEdg, aPT1, aLoc1);
2131     const Handle(Poly_PolygonOnTriangulation)& aPTEF2 = 
2132                                 BRep_Tool::PolygonOnTriangulation(anEdg, aPT2, aLoc2);
2133
2134     TShort_Array1OfShortReal& aNormArr1 = aPT1->ChangeNormals();
2135     TShort_Array1OfShortReal& aNormArr2 = aPT2->ChangeNormals();
2136
2137     if (aPTEF1->Nodes().Lower() != aPTEF2->Nodes().Lower() || 
2138         aPTEF1->Nodes().Upper() != aPTEF2->Nodes().Upper()) 
2139       continue; 
2140
2141     for(Standard_Integer anEdgNode = aPTEF1->Nodes().Lower();
2142                          anEdgNode <= aPTEF1->Nodes().Upper(); anEdgNode++)
2143     {
2144       //Number of node
2145       const Standard_Integer aFNodF1 = aPTEF1->Nodes().Value(anEdgNode);
2146       const Standard_Integer aFNodF2 = aPTEF2->Nodes().Value(anEdgNode);
2147
2148       const Standard_Integer aFNorm1FirstIndex = aNormArr1.Lower() + 3*
2149                                                     (aFNodF1 - aPT1->Nodes().Lower());
2150       const Standard_Integer aFNorm2FirstIndex = aNormArr2.Lower() + 3*
2151                                                     (aFNodF2 - aPT2->Nodes().Lower());
2152
2153       gp_XYZ aNorm1(aNormArr1.Value(aFNorm1FirstIndex),
2154                     aNormArr1.Value(aFNorm1FirstIndex+1),
2155                     aNormArr1.Value(aFNorm1FirstIndex+2));
2156       gp_XYZ aNorm2(aNormArr2.Value(aFNorm2FirstIndex),
2157                     aNormArr2.Value(aFNorm2FirstIndex+1),
2158                     aNormArr2.Value(aFNorm2FirstIndex+2));
2159       const Standard_Real aDot = aNorm1 * aNorm2;
2160
2161       if(aDot > aThresDot)
2162       {
2163         gp_XYZ aNewNorm = (aNorm1 + aNorm2).Normalized();
2164         aNormArr1.ChangeValue(aFNorm1FirstIndex) =
2165           aNormArr2.ChangeValue(aFNorm2FirstIndex) =
2166           static_cast<Standard_ShortReal>(aNewNorm.X());
2167         aNormArr1.ChangeValue(aFNorm1FirstIndex+1) =
2168           aNormArr2.ChangeValue(aFNorm2FirstIndex+1) =
2169           static_cast<Standard_ShortReal>(aNewNorm.Y());
2170         aNormArr1.ChangeValue(aFNorm1FirstIndex+2) =
2171           aNormArr2.ChangeValue(aFNorm2FirstIndex+2) =
2172           static_cast<Standard_ShortReal>(aNewNorm.Z());
2173         aRetVal = Standard_True;
2174       }
2175     }
2176   }
2177
2178   return aRetVal;
2179 }
2180
2181 //=======================================================================
2182 //function : SortFaces
2183 //purpose  : 
2184 //=======================================================================
2185
2186 void  BRepLib::SortFaces (const TopoDS_Shape& Sh,
2187   TopTools_ListOfShape& LF)
2188 {
2189   LF.Clear();
2190   TopTools_ListOfShape LTri,LPlan,LCyl,LCon,LSphere,LTor,LOther;
2191   TopExp_Explorer exp(Sh,TopAbs_FACE);
2192   TopLoc_Location l;
2193   Handle(Geom_Surface) S;
2194
2195   for (; exp.More(); exp.Next()) {
2196     const TopoDS_Face&   F = TopoDS::Face(exp.Current());
2197     S = BRep_Tool::Surface(F, l);
2198     if (!S.IsNull()) {
2199       if (S->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
2200         S = Handle(Geom_RectangularTrimmedSurface)::DownCast (S)->BasisSurface();
2201       }
2202       GeomAdaptor_Surface AS(S);
2203       switch (AS.GetType()) {
2204       case GeomAbs_Plane: 
2205         {
2206           LPlan.Append(F);
2207           break;
2208         }
2209       case GeomAbs_Cylinder: 
2210         {
2211           LCyl.Append(F);
2212           break;
2213         }
2214       case GeomAbs_Cone: 
2215         {
2216           LCon.Append(F);
2217           break;
2218         }
2219       case GeomAbs_Sphere: 
2220         {
2221           LSphere.Append(F);
2222           break;
2223         }
2224       case GeomAbs_Torus: 
2225         {
2226           LTor.Append(F);
2227           break;
2228         }
2229       default:
2230         LOther.Append(F);
2231       }
2232     }
2233     else LTri.Append(F);
2234   }
2235   LF.Append(LPlan); LF.Append(LCyl  ); LF.Append(LCon); LF.Append(LSphere);
2236   LF.Append(LTor ); LF.Append(LOther); LF.Append(LTri); 
2237 }
2238
2239 //=======================================================================
2240 //function : ReverseSortFaces
2241 //purpose  : 
2242 //=======================================================================
2243
2244 void  BRepLib::ReverseSortFaces (const TopoDS_Shape& Sh,
2245   TopTools_ListOfShape& LF)
2246 {
2247   LF.Clear();
2248   // Use the allocator of the result LF for intermediate results
2249   TopTools_ListOfShape LTri(LF.Allocator()), LPlan(LF.Allocator()),
2250     LCyl(LF.Allocator()), LCon(LF.Allocator()), LSphere(LF.Allocator()),
2251     LTor(LF.Allocator()), LOther(LF.Allocator());
2252   TopExp_Explorer exp(Sh,TopAbs_FACE);
2253   TopLoc_Location l;
2254
2255   for (; exp.More(); exp.Next()) {
2256     const TopoDS_Face&   F = TopoDS::Face(exp.Current());
2257     const Handle(Geom_Surface)& S = BRep_Tool::Surface(F, l);
2258     if (!S.IsNull()) {
2259       GeomAdaptor_Surface AS(S);
2260       switch (AS.GetType()) {
2261       case GeomAbs_Plane: 
2262         {
2263           LPlan.Append(F);
2264           break;
2265         }
2266       case GeomAbs_Cylinder: 
2267         {
2268           LCyl.Append(F);
2269           break;
2270         }
2271       case GeomAbs_Cone: 
2272         {
2273           LCon.Append(F);
2274           break;
2275         }
2276       case GeomAbs_Sphere: 
2277         {
2278           LSphere.Append(F);
2279           break;
2280         }
2281       case GeomAbs_Torus: 
2282         {
2283           LTor.Append(F);
2284           break;
2285         }
2286       default:
2287         LOther.Append(F);
2288       }
2289     }
2290     else LTri.Append(F);
2291   }
2292   LF.Append(LTri); LF.Append(LOther); LF.Append(LTor ); LF.Append(LSphere);
2293   LF.Append(LCon); LF.Append(LCyl  ); LF.Append(LPlan);
2294
2295 }
2296
2297 //=======================================================================
2298 // function: BoundingVertex
2299 // purpose : 
2300 //=======================================================================
2301 void BRepLib::BoundingVertex(const NCollection_List<TopoDS_Shape>& theLV,
2302                              gp_Pnt& theNewCenter, Standard_Real& theNewTol)
2303 {
2304   Standard_Integer aNb;
2305   //
2306   aNb=theLV.Extent();
2307   if (aNb < 2) {
2308     return;
2309   }
2310   //
2311   else if (aNb==2) {
2312     Standard_Integer m, n;
2313     Standard_Real aR[2], dR, aD, aEps;
2314     TopoDS_Vertex aV[2];
2315     gp_Pnt aP[2];
2316     //
2317     aEps=RealEpsilon();
2318     for (m=0; m<aNb; ++m) {
2319       aV[m]=(!m)? 
2320         *((TopoDS_Vertex*)(&theLV.First())):
2321         *((TopoDS_Vertex*)(&theLV.Last()));
2322       aP[m]=BRep_Tool::Pnt(aV[m]);
2323       aR[m]=BRep_Tool::Tolerance(aV[m]);
2324     }  
2325     //
2326     m=0; // max R
2327     n=1; // min R
2328     if (aR[0]<aR[1]) {
2329       m=1;
2330       n=0;
2331     }
2332     //
2333     dR=aR[m]-aR[n]; // dR >= 0.
2334     gp_Vec aVD(aP[m], aP[n]);
2335     aD=aVD.Magnitude();
2336     //
2337     if (aD<=dR || aD<aEps) { 
2338       theNewCenter = aP[m];
2339       theNewTol = aR[m];
2340     }
2341     else {
2342       Standard_Real aRr;
2343       gp_XYZ aXYZr;
2344       gp_Pnt aPr;
2345       //
2346       aRr=0.5*(aR[m]+aR[n]+aD);
2347       aXYZr=0.5*(aP[m].XYZ()+aP[n].XYZ()-aVD.XYZ()*(dR/aD));
2348       aPr.SetXYZ(aXYZr);
2349       //
2350       theNewCenter = aPr;
2351       theNewTol = aRr;
2352       //aBB.MakeVertex (aVnew, aPr, aRr);
2353     }
2354     return;
2355   }// else if (aNb==2) {
2356   //
2357   else { // if (aNb>2)
2358     // compute the point
2359     //
2360     // issue 0027540 - sum of doubles may depend on the order
2361     // of addition, thus sort the coordinates for stable result
2362     Standard_Integer i;
2363     NCollection_Array1<gp_Pnt> aPoints(0, aNb-1);
2364     NCollection_List<TopoDS_Edge>::Iterator aIt(theLV);
2365     for (i = 0; aIt.More(); aIt.Next(), ++i) {
2366       const TopoDS_Vertex& aVi = *((TopoDS_Vertex*)(&aIt.Value()));
2367       gp_Pnt aPi = BRep_Tool::Pnt(aVi);
2368       aPoints(i) = aPi;
2369     }
2370     //
2371     std::sort(aPoints.begin(), aPoints.end(), BRepLib_ComparePoints());
2372     //
2373     gp_XYZ aXYZ(0., 0., 0.);
2374     for (i = 0; i < aNb; ++i) {
2375       aXYZ += aPoints(i).XYZ();
2376     }
2377     aXYZ.Divide((Standard_Real)aNb);
2378     //
2379     gp_Pnt aP(aXYZ);
2380     //
2381     // compute the tolerance for the new vertex
2382     Standard_Real aTi, aDi, aDmax;
2383     //
2384     aDmax=-1.;
2385     aIt.Initialize(theLV);
2386     for (; aIt.More(); aIt.Next()) {
2387       TopoDS_Vertex& aVi=*((TopoDS_Vertex*)(&aIt.Value()));
2388       gp_Pnt aPi=BRep_Tool::Pnt(aVi);
2389       aTi=BRep_Tool::Tolerance(aVi);
2390       aDi=aP.SquareDistance(aPi);
2391       aDi=sqrt(aDi);
2392       aDi=aDi+aTi;
2393       if (aDi > aDmax) {
2394         aDmax=aDi;
2395       }
2396     }
2397     //
2398     theNewCenter = aP;
2399     theNewTol = aDmax;
2400   }
2401 }