8e87f259d248329642d0c0feea90062263a43cbd
[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::MapShapesAndUniqueAncestors(aShape, TopAbs_VERTEX, TopAbs_EDGE, parents);
1474   TColStd_MapOfTransient Initialized;
1475   Standard_Integer nbV = parents.Extent();
1476   for (iCur=1; iCur<=nbV; iCur++) {
1477     tol=0;
1478     const TopoDS_Vertex& V = TopoDS::Vertex(parents.FindKey(iCur));
1479     Bnd_Box box;
1480     box.Add(BRep_Tool::Pnt(V));
1481     gp_Pnt p3d;
1482     for (lConx.Initialize(parents(iCur)); lConx.More(); lConx.Next()) {
1483       const TopoDS_Edge& E = TopoDS::Edge(lConx.Value());
1484       tol=Max(tol, BRep_Tool::Tolerance(E));
1485       if(tol > BigTol) continue;
1486       if(!BRep_Tool::SameRange(E)) continue;
1487       Standard_Real par = BRep_Tool::Parameter(V,E);
1488       Handle(BRep_TEdge)& TE = *((Handle(BRep_TEdge)*)&E.TShape());
1489       BRep_ListIteratorOfListOfCurveRepresentation itcr(TE->Curves());
1490       const TopLoc_Location& Eloc = E.Location();
1491       while (itcr.More()) {
1492         // For each CurveRepresentation, check the provided parameter
1493         const Handle(BRep_CurveRepresentation)& cr = itcr.Value();
1494         const TopLoc_Location& loc = cr->Location();
1495         TopLoc_Location L = (Eloc * loc);
1496         if (cr->IsCurve3D()) {
1497           const Handle(Geom_Curve)& C = cr->Curve3D();
1498           if (!C.IsNull()) { // edge non degenerated
1499             p3d = C->Value(par);
1500             p3d.Transform(L.Transformation());
1501             box.Add(p3d);
1502           }
1503         }
1504         else if (cr->IsCurveOnSurface()) {
1505           const Handle(Geom_Surface)& Su = cr->Surface();
1506           const Handle(Geom2d_Curve)& PC = cr->PCurve();
1507           Handle(Geom2d_Curve) PC2;
1508           if (cr->IsCurveOnClosedSurface()) {
1509             PC2 = cr->PCurve2();
1510           }
1511           gp_Pnt2d p2d = PC->Value(par);
1512           p3d = Su->Value(p2d.X(),p2d.Y());
1513           p3d.Transform(L.Transformation());
1514           box.Add(p3d);
1515           if (!PC2.IsNull()) {
1516             p2d = PC2->Value(par);
1517             p3d = Su->Value(p2d.X(),p2d.Y());
1518             p3d.Transform(L.Transformation());
1519             box.Add(p3d);
1520           }
1521         }
1522         itcr.Next();
1523       }
1524     }
1525     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
1526     box.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
1527     aXmax -= aXmin; aYmax -= aYmin; aZmax -= aZmin;
1528     tol = Max(tol,sqrt(aXmax*aXmax+aYmax*aYmax+aZmax*aZmax));
1529     tol += 2.*Epsilon(tol);
1530     if (verifyTolerance) {
1531       // ASet minimum value of the tolerance 
1532       // Attention to sharing of the vertex by other shapes
1533       const Handle(BRep_TVertex)& TV = *((Handle(BRep_TVertex)*)&V.TShape());
1534       if (Initialized.Add(TV)) 
1535         TV->Tolerance(tol);
1536       else 
1537         B.UpdateVertex(V, tol);
1538     }
1539     else {
1540       // Update can only increase tolerance, so if the edge has a greater
1541       //  tolerance than its faces it is not concerned
1542       B.UpdateVertex(V, tol);
1543     }
1544   }
1545 }
1546
1547 //=======================================================================
1548 //function : UpdateInnerTolerances
1549 //purpose  : 
1550 //=======================================================================
1551 void  BRepLib::UpdateInnerTolerances(const TopoDS_Shape& aShape)
1552 {
1553   TopTools_IndexedDataMapOfShapeListOfShape EFmap;
1554   TopExp::MapShapesAndAncestors(aShape, TopAbs_EDGE, TopAbs_FACE, EFmap);
1555   BRep_Builder BB;
1556   for (Standard_Integer i = 1; i <= EFmap.Extent(); i++)
1557   {
1558     TopoDS_Edge anEdge = TopoDS::Edge(EFmap.FindKey(i));
1559     TopoDS_Vertex V1, V2;
1560     TopExp::Vertices(anEdge, V1, V2);
1561     Standard_Real fpar, lpar;
1562     BRep_Tool::Range(anEdge, fpar, lpar);
1563     Standard_Real TolEdge = BRep_Tool::Tolerance(anEdge);
1564     gp_Pnt Pnt1, Pnt2;
1565     Handle(BRepAdaptor_HCurve) anHCurve = new BRepAdaptor_HCurve();
1566     anHCurve->ChangeCurve().Initialize(anEdge);
1567     if (!V1.IsNull())
1568       Pnt1 = BRep_Tool::Pnt(V1);
1569     if (!V2.IsNull())
1570       Pnt2 = BRep_Tool::Pnt(V2);
1571     
1572     if (!BRep_Tool::Degenerated(anEdge) &&
1573         EFmap(i).Extent() > 0)
1574     {
1575       NCollection_Sequence<Handle(Adaptor3d_HCurve)> theRep;
1576       theRep.Append(anHCurve);
1577       TopTools_ListIteratorOfListOfShape itl(EFmap(i));
1578       for (; itl.More(); itl.Next())
1579       {
1580         const TopoDS_Face& aFace = TopoDS::Face(itl.Value());
1581         Handle(BRepAdaptor_HCurve) anHCurvOnSurf = new BRepAdaptor_HCurve();
1582         anHCurvOnSurf->ChangeCurve().Initialize(anEdge, aFace);
1583         theRep.Append(anHCurvOnSurf);
1584       }
1585       
1586       const Standard_Integer NbSamples = (BRep_Tool::SameParameter(anEdge))? 23 : 2;
1587       Standard_Real delta = (lpar - fpar)/(NbSamples-1);
1588       Standard_Real MaxDist = 0.;
1589       for (Standard_Integer j = 2; j <= theRep.Length(); j++)
1590       {
1591         for (Standard_Integer k = 0; k <= NbSamples; k++)
1592         {
1593           Standard_Real ParamOnCenter = (k == NbSamples)? lpar :
1594             fpar + k*delta;
1595           gp_Pnt Center = theRep(1)->Value(ParamOnCenter);
1596           Standard_Real ParamOnCurve = (BRep_Tool::SameParameter(anEdge))? ParamOnCenter
1597             : ((k == 0)? theRep(j)->FirstParameter() : theRep(j)->LastParameter());
1598           gp_Pnt aPoint = theRep(j)->Value(ParamOnCurve);
1599           Standard_Real aDist = Center.Distance(aPoint);
1600           //aDist *= 1.1;
1601           aDist += 2.*Epsilon(aDist);
1602           if (aDist > MaxDist)
1603             MaxDist = aDist;
1604
1605           //Update tolerances of vertices
1606           if (k == 0 && !V1.IsNull())
1607           {
1608             Standard_Real aDist1 = Pnt1.Distance(aPoint);
1609             aDist1 += 2.*Epsilon(aDist1);
1610             BB.UpdateVertex(V1, aDist1);
1611           }
1612           if (k == NbSamples && !V2.IsNull())
1613           {
1614             Standard_Real aDist2 = Pnt2.Distance(aPoint);
1615             aDist2 += 2.*Epsilon(aDist2);
1616             BB.UpdateVertex(V2, aDist2);
1617           }
1618         }
1619       }
1620       BB.UpdateEdge(anEdge, MaxDist);
1621     }
1622     TolEdge = BRep_Tool::Tolerance(anEdge);
1623     if (!V1.IsNull())
1624     {
1625       gp_Pnt End1 = anHCurve->Value(fpar);
1626       Standard_Real dist1 = Pnt1.Distance(End1);
1627       dist1 += 2.*Epsilon(dist1);
1628       BB.UpdateVertex(V1, Max(dist1, TolEdge));
1629     }
1630     if (!V2.IsNull())
1631     {
1632       gp_Pnt End2 = anHCurve->Value(lpar);
1633       Standard_Real dist2 = Pnt2.Distance(End2);
1634       dist2 += 2.*Epsilon(dist2);
1635       BB.UpdateVertex(V2, Max(dist2, TolEdge));
1636     }
1637   }
1638 }
1639
1640 //=======================================================================
1641 //function : OrientClosedSolid
1642 //purpose  : 
1643 //=======================================================================
1644 Standard_Boolean BRepLib::OrientClosedSolid(TopoDS_Solid& solid) 
1645 {
1646   // Set material inside the solid
1647   BRepClass3d_SolidClassifier where(solid);
1648   where.PerformInfinitePoint(Precision::Confusion());
1649   if (where.State()==TopAbs_IN) {
1650     solid.Reverse();
1651   }
1652   else if (where.State()==TopAbs_ON || where.State()==TopAbs_UNKNOWN) 
1653     return Standard_False;
1654
1655   return Standard_True;
1656 }
1657
1658 // Structure for calculation of properties, necessary for decision about continuity
1659 class SurfaceProperties
1660 {
1661 public:
1662   SurfaceProperties(const Handle(Geom_Surface)& theSurface,
1663                     const gp_Trsf&              theSurfaceTrsf,
1664                     const Handle(Geom2d_Curve)& theCurve2D,
1665                     const Standard_Boolean      theReversed)
1666     : mySurfaceProps(theSurface, 2, Precision::Confusion()),
1667       mySurfaceTrsf(theSurfaceTrsf),
1668       myCurve2d(theCurve2D),
1669       myIsReversed(theReversed)
1670   {}
1671
1672   // Calculate derivatives on surface related to the point on curve
1673   void Calculate(const Standard_Real theParamOnCurve)
1674   {
1675     gp_Pnt2d aUV;
1676     myCurve2d->D1(theParamOnCurve, aUV, myCurveTangent);
1677     mySurfaceProps.SetParameters(aUV.X(), aUV.Y());
1678   }
1679
1680   // Returns point just calculated
1681   gp_Pnt Value() 
1682   { return mySurfaceProps.Value().Transformed(mySurfaceTrsf); }
1683
1684   // Calculate a derivative orthogonal to curve's tangent vector
1685   gp_Vec Derivative()
1686   {
1687     gp_Vec aDeriv;
1688     // direction orthogonal to tangent vector of the curve
1689     gp_Vec2d anOrtho(-myCurveTangent.Y(), myCurveTangent.X());
1690     Standard_Real aLen = anOrtho.Magnitude();
1691     if (aLen < Precision::Confusion())
1692       return aDeriv;
1693     anOrtho /= aLen;
1694     if (myIsReversed)
1695       anOrtho.Reverse();
1696
1697     aDeriv.SetLinearForm(anOrtho.X(), mySurfaceProps.D1U(),
1698                          anOrtho.Y(), mySurfaceProps.D1V());
1699     return aDeriv.Transformed(mySurfaceTrsf);
1700   }
1701
1702   // Calculate principal curvatures, which consist of minimal and maximal normal curvatures and
1703   // the directions on the tangent plane (principal direction) where the extremums are reached
1704   void Curvature(gp_Dir& thePrincipalDir1, Standard_Real& theCurvature1,
1705                  gp_Dir& thePrincipalDir2, Standard_Real& theCurvature2)
1706   {
1707     mySurfaceProps.CurvatureDirections(thePrincipalDir1, thePrincipalDir2);
1708     theCurvature1 = mySurfaceProps.MaxCurvature();
1709     theCurvature2 = mySurfaceProps.MinCurvature();
1710     if (myIsReversed)
1711     {
1712       theCurvature1 = -theCurvature1;
1713       theCurvature2 = -theCurvature2;
1714     }
1715     thePrincipalDir1.Transform(mySurfaceTrsf);
1716     thePrincipalDir2.Transform(mySurfaceTrsf);
1717   }
1718
1719 private:
1720   GeomLProp_SLProps    mySurfaceProps; // properties calculator
1721   gp_Trsf              mySurfaceTrsf;
1722   Handle(Geom2d_Curve) myCurve2d;
1723   Standard_Boolean     myIsReversed; // the face based on the surface is reversed
1724
1725   // tangent vector to Pcurve in UV
1726   gp_Vec2d myCurveTangent;
1727 };
1728
1729 //=======================================================================
1730 //function : tgtfaces
1731 //purpose  : check the angle at the border between two squares.
1732 //           Two shares should have a shared front edge.
1733 //=======================================================================
1734 static GeomAbs_Shape tgtfaces(const TopoDS_Edge& Ed,
1735                               const TopoDS_Face& F1,
1736                               const TopoDS_Face& F2,
1737                               const Standard_Real theAngleTol)
1738 {
1739   Standard_Boolean isSeam = F1.IsEqual(F2);
1740
1741   TopoDS_Edge E = Ed;
1742
1743   // Check if pcurves exist on both faces of edge
1744   Standard_Real aFirst,aLast;
1745   E.Orientation(TopAbs_FORWARD);
1746   Handle(Geom2d_Curve) aCurve1 = BRep_Tool::CurveOnSurface(E, F1, aFirst, aLast);
1747   if(aCurve1.IsNull())
1748     return GeomAbs_C0;
1749
1750   if (isSeam)
1751     E.Orientation(TopAbs_REVERSED);
1752   Handle(Geom2d_Curve) aCurve2 = BRep_Tool::CurveOnSurface(E, F2, aFirst, aLast);
1753   if(aCurve2.IsNull())
1754     return GeomAbs_C0;
1755
1756   TopLoc_Location aLoc1, aLoc2;
1757   Handle(Geom_Surface) aSurface1 = BRep_Tool::Surface(F1, aLoc1);
1758   const gp_Trsf& aSurf1Trsf = aLoc1.Transformation();
1759   Handle(Geom_Surface) aSurface2 = BRep_Tool::Surface(F2, aLoc2);
1760   const gp_Trsf& aSurf2Trsf = aLoc2.Transformation();
1761
1762   if (aSurface1->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
1763     aSurface1 = Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface1)->BasisSurface();
1764   if (aSurface2->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
1765     aSurface2 = Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface2)->BasisSurface();
1766
1767   // seam edge on elementary surface is always CN
1768   Standard_Boolean isElementary =
1769     (aSurface1->IsKind(STANDARD_TYPE(Geom_ElementarySurface)) &&
1770      aSurface2->IsKind(STANDARD_TYPE(Geom_ElementarySurface)));
1771   if (isSeam && isElementary)
1772   {
1773     return GeomAbs_CN;
1774   }
1775
1776   SurfaceProperties aSP1(aSurface1, aSurf1Trsf, aCurve1, F1.Orientation() == TopAbs_REVERSED);
1777   SurfaceProperties aSP2(aSurface2, aSurf2Trsf, aCurve2, F2.Orientation() == TopAbs_REVERSED);
1778
1779   Standard_Real f, l, eps;
1780   BRep_Tool::Range(E,f,l);
1781   Extrema_LocateExtPC ext;
1782   Handle(BRepAdaptor_HCurve) aHC2;
1783
1784   eps = (l - f)/100.;
1785   f += eps; // to avoid calculations on  
1786   l -= eps; // points of pointed squares.
1787
1788   const Standard_Real anAngleTol2 = theAngleTol * theAngleTol;
1789
1790   gp_Vec aDer1, aDer2;
1791   gp_Vec aNorm1;
1792   Standard_Real aSqLen1, aSqLen2;
1793   gp_Dir aCrvDir1[2], aCrvDir2[2];
1794   Standard_Real aCrvLen1[2], aCrvLen2[2];
1795
1796   GeomAbs_Shape aCont = (isElementary ? GeomAbs_CN : GeomAbs_C2);
1797   GeomAbs_Shape aCurCont;
1798   Standard_Real u;
1799   for (Standard_Integer i = 0; i <= 20 && aCont > GeomAbs_C0; i++)
1800   {
1801     // First suppose that this is sameParameter
1802     u = f + (l-f)*i/20;
1803
1804     // Check conditions for G1 and C1 continuity:
1805     // * calculate a derivative in tangent plane of each surface
1806     //   orthogonal to curve's tangent vector
1807     // * continuity is C1 if the vectors are equal
1808     // * continuity is G1 if the vectors are just parallel
1809     aCurCont = GeomAbs_C0;
1810
1811     aSP1.Calculate(u);
1812     aSP2.Calculate(u);
1813
1814     aDer1 = aSP1.Derivative();
1815     aSqLen1 = aDer1.SquareMagnitude();
1816     aDer2 = aSP2.Derivative();
1817     aSqLen2 = aDer2.SquareMagnitude();
1818     Standard_Boolean isSmoothSuspect = (aDer1.CrossSquareMagnitude(aDer2) <= anAngleTol2 * aSqLen1 * aSqLen2);
1819     if (!isSmoothSuspect)
1820     {
1821       // Refine by projection
1822       if (aHC2.IsNull())
1823       {
1824         // adaptor for pcurve on the second surface
1825         aHC2 = new BRepAdaptor_HCurve(BRepAdaptor_Curve(E, F2));
1826         ext.Initialize(aHC2->Curve(), f, l, Precision::PConfusion());
1827       }
1828       ext.Perform(aSP1.Value(), u);
1829       if (ext.IsDone() && ext.IsMin())
1830       {
1831         const Extrema_POnCurv& poc = ext.Point();
1832         aSP2.Calculate(poc.Parameter());
1833         aDer2 = aSP2.Derivative();
1834         aSqLen2 = aDer2.SquareMagnitude();
1835       }
1836       isSmoothSuspect = (aDer1.CrossSquareMagnitude(aDer2) <= anAngleTol2 * aSqLen1 * aSqLen2);
1837     }
1838     if (isSmoothSuspect)
1839     {
1840       aCurCont = GeomAbs_G1;
1841       if (Abs(Sqrt(aSqLen1) - Sqrt(aSqLen2)) < Precision::Confusion() &&
1842           aDer1.Dot(aDer2) > Precision::SquareConfusion()) // <= check vectors are codirectional
1843         aCurCont = GeomAbs_C1;
1844     }
1845     else
1846       return GeomAbs_C0;
1847
1848     if (aCont < GeomAbs_G2)
1849       continue; // no need further processing, because maximal continuity is less than G2
1850
1851     // Check conditions for G2 and C2 continuity:
1852     // * calculate principal curvatures on each surface
1853     // * continuity is C2 if directions of principal curvatures are equal on differenct surfaces
1854     // * continuity is G2 if directions of principal curvatures are just parallel
1855     //   and values of curvatures are the same
1856     aSP1.Curvature(aCrvDir1[0], aCrvLen1[0], aCrvDir1[1], aCrvLen1[1]);
1857     aSP2.Curvature(aCrvDir2[0], aCrvLen2[0], aCrvDir2[1], aCrvLen2[1]);
1858     for (Standard_Integer aStep = 0; aStep <= 1; ++aStep)
1859     {
1860       if (aCrvDir1[0].XYZ().CrossSquareMagnitude(aCrvDir2[aStep].XYZ()) <= Precision::SquareConfusion() &&
1861           Abs(aCrvLen1[0] - aCrvLen2[aStep]) < Precision::Confusion() &&
1862           aCrvDir1[1].XYZ().CrossSquareMagnitude(aCrvDir2[1 - aStep].XYZ()) <= Precision::SquareConfusion() &&
1863           Abs(aCrvLen1[1] - aCrvLen2[1 - aStep]) < Precision::Confusion())
1864       {
1865         if (aCurCont == GeomAbs_C1 &&
1866             aCrvDir1[0].Dot(aCrvDir2[aStep]) > Precision::Confusion() &&
1867             aCrvDir1[1].Dot(aCrvDir2[1 - aStep]) > Precision::Confusion())
1868           aCurCont = GeomAbs_C2;
1869         else
1870           aCurCont = GeomAbs_G2;
1871         break;
1872       }
1873     }
1874
1875     if (aCurCont < aCont)
1876       aCont = aCurCont;
1877   }
1878
1879   // according to the list of supported elementary surfaces,
1880   // if the continuity is C2, than it is totally CN
1881   if (isElementary && aCont == GeomAbs_C2)
1882     aCont = GeomAbs_CN;
1883   return aCont;
1884 }
1885
1886 //=======================================================================
1887 // function : EncodeRegularity
1888 // purpose  : Code the regularities on all edges of the shape, boundary of 
1889 //            two faces that do not have it.
1890 //            Takes into account that compound may consists of same solid
1891 //            placed with different transformations
1892 //=======================================================================
1893 static void EncodeRegularity(const TopoDS_Shape&        theShape,
1894                              const Standard_Real        theTolAng,
1895                              TopTools_MapOfShape&       theMap,
1896                              const TopTools_MapOfShape& theEdgesToEncode = TopTools_MapOfShape())
1897 {
1898   TopoDS_Shape aShape = theShape;
1899   TopLoc_Location aNullLoc;
1900   aShape.Location(aNullLoc); // nullify location
1901   if (!theMap.Add(aShape))
1902     return; // do not need to process shape twice
1903
1904   if (aShape.ShapeType() == TopAbs_COMPOUND ||
1905       aShape.ShapeType() == TopAbs_COMPSOLID)
1906   {
1907     for (TopoDS_Iterator it(aShape); it.More(); it.Next())
1908       EncodeRegularity(it.Value(), theTolAng, theMap, theEdgesToEncode);
1909     return;
1910   }
1911
1912   try {
1913     OCC_CATCH_SIGNALS
1914
1915     TopTools_IndexedDataMapOfShapeListOfShape M;
1916     TopExp::MapShapesAndAncestors(aShape, TopAbs_EDGE, TopAbs_FACE, M);
1917     TopTools_ListIteratorOfListOfShape It;
1918     TopExp_Explorer Ex;
1919     TopoDS_Face F1,F2;
1920     Standard_Boolean found;
1921     for (Standard_Integer i = 1; i <= M.Extent(); i++){
1922       TopoDS_Edge E = TopoDS::Edge(M.FindKey(i));
1923       if (!theEdgesToEncode.IsEmpty())
1924       {
1925         // process only the edges from the list to update their regularity
1926         TopoDS_Shape aPureEdge = E.Located(aNullLoc);
1927         aPureEdge.Orientation(TopAbs_FORWARD);
1928         if (!theEdgesToEncode.Contains(aPureEdge))
1929           continue;
1930       }
1931
1932       found = Standard_False;                                     
1933       F1.Nullify();
1934       for (It.Initialize(M.FindFromIndex(i)); It.More() && !found; It.Next()){
1935         if (F1.IsNull()) { F1 = TopoDS::Face(It.Value()); }
1936         else {
1937           const TopoDS_Face& aTmpF2 = TopoDS::Face(It.Value());
1938           if (!F1.IsSame(aTmpF2)){
1939             found = Standard_True;
1940             F2 = aTmpF2;
1941           }
1942         }
1943       }
1944       if (!found && !F1.IsNull()){//is it a sewing edge?
1945         TopAbs_Orientation orE = E.Orientation();
1946         TopoDS_Edge curE;
1947         for (Ex.Init(F1, TopAbs_EDGE); Ex.More() && !found; Ex.Next()){
1948           curE = TopoDS::Edge(Ex.Current());
1949           if (E.IsSame(curE) && orE != curE.Orientation()) {
1950             found = Standard_True;
1951             F2 = F1;
1952           }
1953         }
1954       }
1955       if (found)
1956         BRepLib::EncodeRegularity(E, F1, F2, theTolAng);
1957     }
1958   }
1959   catch (Standard_Failure const& anException) {
1960 #ifdef OCCT_DEBUG
1961     cout << "Warning: Exception in BRepLib::EncodeRegularity(): ";
1962     anException.Print(cout);
1963     cout << endl;
1964 #endif
1965     (void)anException;
1966   }
1967 }
1968
1969 //=======================================================================
1970 // function : EncodeRegularity
1971 // purpose  : code the regularities on all edges of the shape, boundary of 
1972 //            two faces that do not have it.
1973 //=======================================================================
1974
1975 void BRepLib::EncodeRegularity(const TopoDS_Shape& S,
1976   const Standard_Real TolAng)
1977 {
1978   TopTools_MapOfShape aMap;
1979   ::EncodeRegularity(S, TolAng, aMap);
1980 }
1981
1982 //=======================================================================
1983 // function : EncodeRegularity
1984 // purpose  : code the regularities on all edges in the list that do not 
1985 //            have it, and which are boundary of two faces on the shape.
1986 //=======================================================================
1987
1988 void BRepLib::EncodeRegularity(const TopoDS_Shape& S,
1989   const TopTools_ListOfShape& LE,
1990   const Standard_Real TolAng)
1991 {
1992   // Collect edges without location and orientation
1993   TopTools_MapOfShape aPureEdges;
1994   TopLoc_Location aNullLoc;
1995   TopTools_ListIteratorOfListOfShape anEdgeIt(LE);
1996   for (; anEdgeIt.More(); anEdgeIt.Next())
1997   {
1998     TopoDS_Shape anEdge = anEdgeIt.Value();
1999     anEdge.Location(aNullLoc);
2000     anEdge.Orientation(TopAbs_FORWARD);
2001     aPureEdges.Add(anEdge);
2002   }
2003
2004   TopTools_MapOfShape aMap;
2005   ::EncodeRegularity(S, TolAng, aMap, aPureEdges);
2006 }
2007
2008 //=======================================================================
2009 // function : EncodeRegularity
2010 // purpose  : code the regularity between 2 faces connected by edge 
2011 //=======================================================================
2012
2013 void BRepLib::EncodeRegularity(TopoDS_Edge& E,
2014   const TopoDS_Face& F1,
2015   const TopoDS_Face& F2,
2016   const Standard_Real TolAng)
2017 {
2018   BRep_Builder B;
2019   if(BRep_Tool::Continuity(E,F1,F2)<=GeomAbs_C0){
2020     try {
2021       GeomAbs_Shape aCont = tgtfaces(E, F1, F2, TolAng);
2022       B.Continuity(E,F1,F2,aCont);
2023       
2024     }
2025     catch(Standard_Failure)
2026     {
2027 #ifdef OCCT_DEBUG
2028       cout << "Failure: Exception in BRepLib::EncodeRegularity" << endl;
2029 #endif
2030     }
2031   }
2032 }
2033
2034 //=======================================================================
2035 // function : EnsureNormalConsistency
2036 // purpose  : Corrects the normals in Poly_Triangulation of faces.
2037 //              Returns TRUE if any correction is done.
2038 //=======================================================================
2039 Standard_Boolean BRepLib::
2040             EnsureNormalConsistency(const TopoDS_Shape& theShape,
2041                                     const Standard_Real theAngTol,
2042                                     const Standard_Boolean theForceComputeNormals)
2043 {
2044   const Standard_Real aThresDot = cos(theAngTol);
2045
2046   Standard_Boolean aRetVal = Standard_False, isNormalsFound = Standard_False;
2047
2048   // compute normals if they are absent
2049   TopExp_Explorer anExpFace(theShape,TopAbs_FACE);
2050   for (; anExpFace.More(); anExpFace.Next())
2051   {
2052     const TopoDS_Face& aFace = TopoDS::Face(anExpFace.Current());
2053     const Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
2054     if(aSurf.IsNull())
2055       continue;
2056     TopLoc_Location aLoc;
2057     const Handle(Poly_Triangulation)& aPT = BRep_Tool::Triangulation(aFace, aLoc);
2058     if(aPT.IsNull())
2059       continue;
2060     if (!theForceComputeNormals && aPT->HasNormals())
2061     {
2062       isNormalsFound = Standard_True;
2063       continue;
2064     }
2065
2066     GeomLProp_SLProps aSLP(aSurf, 2, Precision::Confusion());
2067     const Standard_Integer anArrDim = 3*aPT->NbNodes();
2068     Handle(TShort_HArray1OfShortReal) aNormArr = new TShort_HArray1OfShortReal(1, anArrDim);
2069     Standard_Integer anNormInd = aNormArr->Lower();
2070     for(Standard_Integer i = aPT->UVNodes().Lower(); i <= aPT->UVNodes().Upper(); i++)
2071     {
2072       const gp_Pnt2d &aP2d = aPT->UVNodes().Value(i);
2073       aSLP.SetParameters(aP2d.X(), aP2d.Y());
2074
2075       gp_XYZ aNorm(0.,0.,0.);
2076       if(!aSLP.IsNormalDefined())
2077       {
2078 #ifdef OCCT_DEBUG
2079         cout << "BRepLib::EnsureNormalConsistency(): Cannot find normal!" << endl;
2080 #endif
2081       }
2082       else
2083       {
2084         aNorm = aSLP.Normal().XYZ();
2085         if (aFace.Orientation() == TopAbs_REVERSED)
2086           aNorm.Reverse();
2087       }
2088       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.X());
2089       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.Y());
2090       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.Z());
2091     }
2092
2093     aRetVal = Standard_True;
2094     isNormalsFound = Standard_True;
2095     aPT->SetNormals(aNormArr);
2096   }
2097
2098   if(!isNormalsFound)
2099   {
2100     return aRetVal;
2101   }
2102
2103   // loop by edges
2104   TopTools_IndexedDataMapOfShapeListOfShape aMapEF;
2105   TopExp::MapShapesAndAncestors(theShape,TopAbs_EDGE,TopAbs_FACE,aMapEF);
2106   for(Standard_Integer anInd = 1; anInd <= aMapEF.Extent(); anInd++)
2107   {
2108     const TopoDS_Edge& anEdg = TopoDS::Edge(aMapEF.FindKey(anInd));
2109     const TopTools_ListOfShape& anEdgList = aMapEF.FindFromIndex(anInd);
2110     if (anEdgList.Extent() != 2)
2111       continue;
2112     TopTools_ListIteratorOfListOfShape anItF(anEdgList);
2113     const TopoDS_Face aFace1 = TopoDS::Face(anItF.Value());
2114     anItF.Next();
2115     const TopoDS_Face aFace2 = TopoDS::Face(anItF.Value());
2116     TopLoc_Location aLoc1, aLoc2;
2117     const Handle(Poly_Triangulation)& aPT1 = BRep_Tool::Triangulation(aFace1, aLoc1);
2118     const Handle(Poly_Triangulation)& aPT2 = BRep_Tool::Triangulation(aFace2, aLoc2);
2119
2120     if(aPT1.IsNull() || aPT2.IsNull())
2121       continue;
2122
2123     if(!aPT1->HasNormals() || !aPT2->HasNormals())
2124       continue;
2125
2126     const Handle(Poly_PolygonOnTriangulation)& aPTEF1 = 
2127                                 BRep_Tool::PolygonOnTriangulation(anEdg, aPT1, aLoc1);
2128     const Handle(Poly_PolygonOnTriangulation)& aPTEF2 = 
2129                                 BRep_Tool::PolygonOnTriangulation(anEdg, aPT2, aLoc2);
2130
2131     TShort_Array1OfShortReal& aNormArr1 = aPT1->ChangeNormals();
2132     TShort_Array1OfShortReal& aNormArr2 = aPT2->ChangeNormals();
2133
2134     if (aPTEF1->Nodes().Lower() != aPTEF2->Nodes().Lower() || 
2135         aPTEF1->Nodes().Upper() != aPTEF2->Nodes().Upper()) 
2136       continue; 
2137
2138     for(Standard_Integer anEdgNode = aPTEF1->Nodes().Lower();
2139                          anEdgNode <= aPTEF1->Nodes().Upper(); anEdgNode++)
2140     {
2141       //Number of node
2142       const Standard_Integer aFNodF1 = aPTEF1->Nodes().Value(anEdgNode);
2143       const Standard_Integer aFNodF2 = aPTEF2->Nodes().Value(anEdgNode);
2144
2145       const Standard_Integer aFNorm1FirstIndex = aNormArr1.Lower() + 3*
2146                                                     (aFNodF1 - aPT1->Nodes().Lower());
2147       const Standard_Integer aFNorm2FirstIndex = aNormArr2.Lower() + 3*
2148                                                     (aFNodF2 - aPT2->Nodes().Lower());
2149
2150       gp_XYZ aNorm1(aNormArr1.Value(aFNorm1FirstIndex),
2151                     aNormArr1.Value(aFNorm1FirstIndex+1),
2152                     aNormArr1.Value(aFNorm1FirstIndex+2));
2153       gp_XYZ aNorm2(aNormArr2.Value(aFNorm2FirstIndex),
2154                     aNormArr2.Value(aFNorm2FirstIndex+1),
2155                     aNormArr2.Value(aFNorm2FirstIndex+2));
2156       const Standard_Real aDot = aNorm1 * aNorm2;
2157
2158       if(aDot > aThresDot)
2159       {
2160         gp_XYZ aNewNorm = (aNorm1 + aNorm2).Normalized();
2161         aNormArr1.ChangeValue(aFNorm1FirstIndex) =
2162           aNormArr2.ChangeValue(aFNorm2FirstIndex) =
2163           static_cast<Standard_ShortReal>(aNewNorm.X());
2164         aNormArr1.ChangeValue(aFNorm1FirstIndex+1) =
2165           aNormArr2.ChangeValue(aFNorm2FirstIndex+1) =
2166           static_cast<Standard_ShortReal>(aNewNorm.Y());
2167         aNormArr1.ChangeValue(aFNorm1FirstIndex+2) =
2168           aNormArr2.ChangeValue(aFNorm2FirstIndex+2) =
2169           static_cast<Standard_ShortReal>(aNewNorm.Z());
2170         aRetVal = Standard_True;
2171       }
2172     }
2173   }
2174
2175   return aRetVal;
2176 }
2177
2178 //=======================================================================
2179 //function : SortFaces
2180 //purpose  : 
2181 //=======================================================================
2182
2183 void  BRepLib::SortFaces (const TopoDS_Shape& Sh,
2184   TopTools_ListOfShape& LF)
2185 {
2186   LF.Clear();
2187   TopTools_ListOfShape LTri,LPlan,LCyl,LCon,LSphere,LTor,LOther;
2188   TopExp_Explorer exp(Sh,TopAbs_FACE);
2189   TopLoc_Location l;
2190   Handle(Geom_Surface) S;
2191
2192   for (; exp.More(); exp.Next()) {
2193     const TopoDS_Face&   F = TopoDS::Face(exp.Current());
2194     S = BRep_Tool::Surface(F, l);
2195     if (!S.IsNull()) {
2196       if (S->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
2197         S = Handle(Geom_RectangularTrimmedSurface)::DownCast (S)->BasisSurface();
2198       }
2199       GeomAdaptor_Surface AS(S);
2200       switch (AS.GetType()) {
2201       case GeomAbs_Plane: 
2202         {
2203           LPlan.Append(F);
2204           break;
2205         }
2206       case GeomAbs_Cylinder: 
2207         {
2208           LCyl.Append(F);
2209           break;
2210         }
2211       case GeomAbs_Cone: 
2212         {
2213           LCon.Append(F);
2214           break;
2215         }
2216       case GeomAbs_Sphere: 
2217         {
2218           LSphere.Append(F);
2219           break;
2220         }
2221       case GeomAbs_Torus: 
2222         {
2223           LTor.Append(F);
2224           break;
2225         }
2226       default:
2227         LOther.Append(F);
2228       }
2229     }
2230     else LTri.Append(F);
2231   }
2232   LF.Append(LPlan); LF.Append(LCyl  ); LF.Append(LCon); LF.Append(LSphere);
2233   LF.Append(LTor ); LF.Append(LOther); LF.Append(LTri); 
2234 }
2235
2236 //=======================================================================
2237 //function : ReverseSortFaces
2238 //purpose  : 
2239 //=======================================================================
2240
2241 void  BRepLib::ReverseSortFaces (const TopoDS_Shape& Sh,
2242   TopTools_ListOfShape& LF)
2243 {
2244   LF.Clear();
2245   // Use the allocator of the result LF for intermediate results
2246   TopTools_ListOfShape LTri(LF.Allocator()), LPlan(LF.Allocator()),
2247     LCyl(LF.Allocator()), LCon(LF.Allocator()), LSphere(LF.Allocator()),
2248     LTor(LF.Allocator()), LOther(LF.Allocator());
2249   TopExp_Explorer exp(Sh,TopAbs_FACE);
2250   TopLoc_Location l;
2251
2252   for (; exp.More(); exp.Next()) {
2253     const TopoDS_Face&   F = TopoDS::Face(exp.Current());
2254     const Handle(Geom_Surface)& S = BRep_Tool::Surface(F, l);
2255     if (!S.IsNull()) {
2256       GeomAdaptor_Surface AS(S);
2257       switch (AS.GetType()) {
2258       case GeomAbs_Plane: 
2259         {
2260           LPlan.Append(F);
2261           break;
2262         }
2263       case GeomAbs_Cylinder: 
2264         {
2265           LCyl.Append(F);
2266           break;
2267         }
2268       case GeomAbs_Cone: 
2269         {
2270           LCon.Append(F);
2271           break;
2272         }
2273       case GeomAbs_Sphere: 
2274         {
2275           LSphere.Append(F);
2276           break;
2277         }
2278       case GeomAbs_Torus: 
2279         {
2280           LTor.Append(F);
2281           break;
2282         }
2283       default:
2284         LOther.Append(F);
2285       }
2286     }
2287     else LTri.Append(F);
2288   }
2289   LF.Append(LTri); LF.Append(LOther); LF.Append(LTor ); LF.Append(LSphere);
2290   LF.Append(LCon); LF.Append(LCyl  ); LF.Append(LPlan);
2291
2292 }
2293
2294 //=======================================================================
2295 // function: BoundingVertex
2296 // purpose : 
2297 //=======================================================================
2298 void BRepLib::BoundingVertex(const NCollection_List<TopoDS_Shape>& theLV,
2299                              gp_Pnt& theNewCenter, Standard_Real& theNewTol)
2300 {
2301   Standard_Integer aNb;
2302   //
2303   aNb=theLV.Extent();
2304   if (aNb < 2) {
2305     return;
2306   }
2307   //
2308   else if (aNb==2) {
2309     Standard_Integer m, n;
2310     Standard_Real aR[2], dR, aD, aEps;
2311     TopoDS_Vertex aV[2];
2312     gp_Pnt aP[2];
2313     //
2314     aEps=RealEpsilon();
2315     for (m=0; m<aNb; ++m) {
2316       aV[m]=(!m)? 
2317         *((TopoDS_Vertex*)(&theLV.First())):
2318         *((TopoDS_Vertex*)(&theLV.Last()));
2319       aP[m]=BRep_Tool::Pnt(aV[m]);
2320       aR[m]=BRep_Tool::Tolerance(aV[m]);
2321     }  
2322     //
2323     m=0; // max R
2324     n=1; // min R
2325     if (aR[0]<aR[1]) {
2326       m=1;
2327       n=0;
2328     }
2329     //
2330     dR=aR[m]-aR[n]; // dR >= 0.
2331     gp_Vec aVD(aP[m], aP[n]);
2332     aD=aVD.Magnitude();
2333     //
2334     if (aD<=dR || aD<aEps) { 
2335       theNewCenter = aP[m];
2336       theNewTol = aR[m];
2337     }
2338     else {
2339       Standard_Real aRr;
2340       gp_XYZ aXYZr;
2341       gp_Pnt aPr;
2342       //
2343       aRr=0.5*(aR[m]+aR[n]+aD);
2344       aXYZr=0.5*(aP[m].XYZ()+aP[n].XYZ()-aVD.XYZ()*(dR/aD));
2345       aPr.SetXYZ(aXYZr);
2346       //
2347       theNewCenter = aPr;
2348       theNewTol = aRr;
2349       //aBB.MakeVertex (aVnew, aPr, aRr);
2350     }
2351     return;
2352   }// else if (aNb==2) {
2353   //
2354   else { // if (aNb>2)
2355     // compute the point
2356     //
2357     // issue 0027540 - sum of doubles may depend on the order
2358     // of addition, thus sort the coordinates for stable result
2359     Standard_Integer i;
2360     NCollection_Array1<gp_Pnt> aPoints(0, aNb-1);
2361     NCollection_List<TopoDS_Edge>::Iterator aIt(theLV);
2362     for (i = 0; aIt.More(); aIt.Next(), ++i) {
2363       const TopoDS_Vertex& aVi = *((TopoDS_Vertex*)(&aIt.Value()));
2364       gp_Pnt aPi = BRep_Tool::Pnt(aVi);
2365       aPoints(i) = aPi;
2366     }
2367     //
2368     std::sort(aPoints.begin(), aPoints.end(), BRepLib_ComparePoints());
2369     //
2370     gp_XYZ aXYZ(0., 0., 0.);
2371     for (i = 0; i < aNb; ++i) {
2372       aXYZ += aPoints(i).XYZ();
2373     }
2374     aXYZ.Divide((Standard_Real)aNb);
2375     //
2376     gp_Pnt aP(aXYZ);
2377     //
2378     // compute the tolerance for the new vertex
2379     Standard_Real aTi, aDi, aDmax;
2380     //
2381     aDmax=-1.;
2382     aIt.Initialize(theLV);
2383     for (; aIt.More(); aIt.Next()) {
2384       TopoDS_Vertex& aVi=*((TopoDS_Vertex*)(&aIt.Value()));
2385       gp_Pnt aPi=BRep_Tool::Pnt(aVi);
2386       aTi=BRep_Tool::Tolerance(aVi);
2387       aDi=aP.SquareDistance(aPi);
2388       aDi=sqrt(aDi);
2389       aDi=aDi+aTi;
2390       if (aDi > aDmax) {
2391         aDmax=aDi;
2392       }
2393     }
2394     //
2395     theNewCenter = aP;
2396     theNewTol = aDmax;
2397   }
2398 }