0028643: Coding rules - eliminate GCC compiler warnings -Wmisleading-indentation
[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   for (Standard_Integer i = 0; i <= nbp; ++i)
920   {
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     dist(i+1) = temp;
962
963     d2 = Max (d2, temp);
964   }
965
966   if(Precision::IsInfinite(d2))
967   {
968     return d2;
969   }
970
971   d2 = Sqrt(d2);
972   if(dapp > d2)
973   {
974     return dapp;
975   }
976
977   Standard_Boolean ana = Standard_False;
978   Standard_Real D2 = 0;
979   Standard_Integer N1 = 0;
980   Standard_Integer N2 = 0;
981   Standard_Integer N3 = 0;
982
983   for (Standard_Integer i = 1; i<= nbp+10; ++i)
984   {
985     if (dist(i) > 0)
986     {
987       if (dist(i) < 1.0)
988       {
989         ++N1;
990       }
991       else
992       {
993         ++N2;
994       }
995     }
996   }
997
998   if (N1 > N2 && N2 != 0)
999   {
1000     N3 = 100*N2/(N1+N2);
1001   }
1002   if (N3 < 10 && N3 != 0)
1003   {
1004     ana = Standard_True;
1005     for (Standard_Integer i = 1; i <= nbp+10; ++i)
1006     {
1007       if (dist(i) > 0 && dist(i) < 1.0)
1008       {
1009         D2 = Max (D2, dist(i));
1010       }
1011     }
1012   }
1013
1014   //d2 = 1.5*sqrt(d2);
1015   d2 = (!ana) ? 1.5 * d2 : 1.5*sqrt(D2);
1016   d2 = Max (d2, 1.e-7);
1017   return d2;
1018 }
1019
1020 void BRepLib::SameParameter(const TopoDS_Edge&  AnEdge,
1021   const Standard_Real Tolerance)
1022 {
1023   if (BRep_Tool::SameParameter(AnEdge)) return;
1024
1025   const Standard_Integer NCONTROL = 22;
1026
1027   Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve();
1028   Handle(Geom2dAdaptor_HCurve) HC2d = new Geom2dAdaptor_HCurve();
1029   Handle(GeomAdaptor_HSurface) HS = new GeomAdaptor_HSurface();
1030   GeomAdaptor_Curve& GAC = HC->ChangeCurve();
1031   Geom2dAdaptor_Curve& GAC2d = HC2d->ChangeCurve2d();
1032   GeomAdaptor_Surface& GAS = HS->ChangeSurface();
1033   Standard_Real f3d =0.,l3d =0.;
1034   TopLoc_Location L3d;
1035   Handle(Geom_Curve) C3d;
1036
1037   const Handle(BRep_TEdge)& TE = *((Handle(BRep_TEdge)*) &AnEdge.TShape());
1038   BRep_ListOfCurveRepresentation& CList = TE->ChangeCurves();
1039   BRep_ListIteratorOfListOfCurveRepresentation It(CList);
1040
1041   Standard_Boolean NotDone = Standard_True;
1042
1043   while (NotDone && It.More()) {
1044     Handle(BRep_GCurve) GCurve = Handle(BRep_GCurve)::DownCast(It.Value());
1045     if (!GCurve.IsNull() && GCurve->IsCurve3D()) {
1046       C3d = GCurve->Curve3D() ;
1047       f3d = GCurve->First();
1048       l3d = GCurve->Last();
1049       L3d = GCurve->Location() ;
1050       NotDone = Standard_False;
1051     } 
1052     It.Next() ;
1053   }
1054
1055   if(C3d.IsNull()) return;
1056
1057   // modified by NIZHNY-OCC486  Tue Aug 27 17:15:13 2002 :
1058   Standard_Boolean m_TrimmedPeriodical = Standard_False;
1059   Handle(Standard_Type) TheType = C3d->DynamicType();
1060   if( TheType == STANDARD_TYPE(Geom_TrimmedCurve))
1061   {
1062     Handle(Geom_Curve) gtC (Handle(Geom_TrimmedCurve)::DownCast (C3d)->BasisCurve());
1063     m_TrimmedPeriodical = gtC->IsPeriodic();
1064   }
1065   // modified by NIZHNY-OCC486  Tue Aug 27 17:15:17 2002 .
1066
1067   BRep_Builder B;
1068   if(!C3d->IsPeriodic()) {
1069     Standard_Real Udeb = C3d->FirstParameter();
1070     Standard_Real Ufin = C3d->LastParameter();
1071     // modified by NIZHNY-OCC486  Tue Aug 27 17:17:14 2002 :
1072     //if (Udeb > f3d) f3d = Udeb;
1073     //if (l3d > Ufin) l3d = Ufin;
1074     if(!m_TrimmedPeriodical)
1075     {
1076       if (Udeb > f3d) f3d = Udeb;
1077       if (l3d > Ufin) l3d = Ufin;
1078     }
1079     // modified by NIZHNY-OCC486  Tue Aug 27 17:17:55 2002 .
1080   }
1081   if(!L3d.IsIdentity()){
1082     C3d = Handle(Geom_Curve)::DownCast(C3d->Transformed(L3d.Transformation()));
1083   }
1084   GAC.Load(C3d,f3d,l3d);
1085
1086   Standard_Boolean IsSameP = 1;
1087   Standard_Real maxdist = 0.;
1088
1089   //  Modified by skv - Thu Jun  3 12:39:19 2004 OCC5898 Begin
1090   Standard_Real anEdgeTol = BRep_Tool::Tolerance(AnEdge);
1091   //  Modified by skv - Thu Jun  3 12:39:20 2004 OCC5898 End
1092   Standard_Boolean SameRange = BRep_Tool::SameRange(AnEdge);
1093   Standard_Boolean YaPCu = Standard_False;
1094   const Standard_Real BigError = 1.e10;
1095   It.Initialize(CList);
1096
1097   while (It.More()) {
1098     Standard_Boolean isANA = Standard_False;
1099     Standard_Boolean isBSP = Standard_False;
1100     Handle(BRep_GCurve) GCurve = Handle(BRep_GCurve)::DownCast(It.Value());
1101     Handle(Geom2d_Curve) PC[2];
1102     Handle(Geom_Surface) S;
1103     if (!GCurve.IsNull() && GCurve->IsCurveOnSurface()) {
1104       YaPCu = Standard_True;
1105       PC[0] = GCurve->PCurve();
1106       TopLoc_Location PCLoc = GCurve->Location();
1107       S = GCurve->Surface();
1108       if (!PCLoc.IsIdentity() ) {
1109         S = Handle(Geom_Surface)::DownCast(S->Transformed(PCLoc.Transformation()));
1110       }
1111
1112       GAS.Load(S);
1113       if (GCurve->IsCurveOnClosedSurface()) {
1114         PC[1] = GCurve->PCurve2();
1115       }
1116
1117       // Eval tol2d to compute SameRange
1118       Standard_Real UResol = Max(GAS.UResolution(Tolerance), Precision::PConfusion());
1119       Standard_Real VResol = Max(GAS.VResolution(Tolerance), Precision::PConfusion());
1120       Standard_Real Tol2d  = Min(UResol, VResol);
1121       for(Standard_Integer i = 0; i < 2; i++){
1122         Handle(Geom2d_Curve) curPC = PC[i];
1123         Standard_Boolean updatepc = 0;
1124         if(curPC.IsNull()) break;
1125         if(!SameRange){
1126           GeomLib::SameRange(Tol2d,
1127             PC[i],GCurve->First(),GCurve->Last(),
1128             f3d,l3d,curPC);
1129
1130           updatepc = (curPC != PC[i]);
1131
1132         }
1133         Standard_Boolean goodpc = 1;
1134         GAC2d.Load(curPC,f3d,l3d);
1135
1136         Standard_Real error = ComputeTol(HC, HC2d, HS, NCONTROL);
1137
1138         if(error > BigError)
1139         {
1140           maxdist = error;
1141           break;
1142         }
1143
1144         if(GAC2d.GetType() == GeomAbs_BSplineCurve && 
1145           GAC2d.Continuity() == GeomAbs_C0) {
1146             Handle(Geom2d_BSplineCurve) bs2d = GAC2d.BSpline();
1147             Handle(Geom2d_BSplineCurve) bs2dsov = bs2d;
1148             Standard_Real fC0 = bs2d->FirstParameter(), lC0 = bs2d->LastParameter();
1149             Standard_Boolean repar = Standard_True;
1150             gp_Pnt2d OriginPoint;
1151             bs2d->D0(fC0, OriginPoint);
1152             Geom2dConvert::C0BSplineToC1BSplineCurve(bs2d, Tol2d);
1153             isBSP = Standard_True; 
1154
1155             if(bs2d->IsPeriodic()) { // -------- IFV, Jan 2000
1156               gp_Pnt2d NewOriginPoint;
1157               bs2d->D0(bs2d->FirstParameter(), NewOriginPoint);
1158               if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1159                 Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) {
1160
1161                   TColStd_Array1OfReal Knotbs2d (1, bs2d->NbKnots());
1162                   bs2d->Knots(Knotbs2d);
1163
1164                   for(Standard_Integer Index = 1; Index <= bs2d->NbKnots(); Index++) {
1165                     bs2d->D0(Knotbs2d(Index), NewOriginPoint);
1166                     if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1167                       Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) continue;
1168
1169                     bs2d->SetOrigin(Index);
1170                     break;
1171                   }
1172               }
1173             }
1174
1175             if(bs2d->Continuity() == GeomAbs_C0) {
1176               Standard_Real tolbail;
1177               if(EvalTol(curPC,S,GAC,Tolerance,tolbail)){
1178                 bs2d = bs2dsov;
1179                 Standard_Real UResbail = GAS.UResolution(tolbail);
1180                 Standard_Real VResbail = GAS.VResolution(tolbail);
1181                 Standard_Real Tol2dbail  = Min(UResbail,VResbail);
1182                 bs2d->D0(bs2d->FirstParameter(), OriginPoint); 
1183
1184                 Standard_Integer nbp = bs2d->NbPoles();
1185                 TColgp_Array1OfPnt2d poles(1,nbp);
1186                 bs2d->Poles(poles);
1187                 gp_Pnt2d p = poles(1), p1;
1188                 Standard_Real d = Precision::Infinite();
1189                 for(Standard_Integer ip = 2; ip <= nbp; ip++) {
1190                   p1 = poles(ip);
1191                   d = Min(d,p.SquareDistance(p1));
1192                   p = p1;
1193                 }
1194                 d = sqrt(d)*.1;
1195
1196                 Tol2dbail = Max(Min(Tol2dbail,d),Tol2d);
1197
1198                 Geom2dConvert::C0BSplineToC1BSplineCurve(bs2d,Tol2dbail);
1199
1200                 if(bs2d->IsPeriodic()) { // -------- IFV, Jan 2000
1201                   gp_Pnt2d NewOriginPoint;
1202                   bs2d->D0(bs2d->FirstParameter(), NewOriginPoint);
1203                   if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1204                     Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) {
1205
1206                       TColStd_Array1OfReal Knotbs2d (1, bs2d->NbKnots());
1207                       bs2d->Knots(Knotbs2d);
1208
1209                       for(Standard_Integer Index = 1; Index <= bs2d->NbKnots(); Index++) {
1210                         bs2d->D0(Knotbs2d(Index), NewOriginPoint);
1211                         if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1212                           Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) continue;
1213
1214                         bs2d->SetOrigin(Index);
1215                         break;
1216                       }
1217                   }
1218                 }
1219
1220
1221                 if(bs2d->Continuity() == GeomAbs_C0) {
1222                   goodpc = 1;
1223                   bs2d = bs2dsov;
1224                   repar = Standard_False;
1225                 }
1226               }
1227               else goodpc = 0;
1228             }
1229
1230             if(goodpc){
1231               if(repar) {
1232                 Standard_Integer NbKnots = bs2d->NbKnots();
1233                 TColStd_Array1OfReal Knots(1,NbKnots);
1234                 bs2d->Knots(Knots);
1235                 //          BSplCLib::Reparametrize(f3d,l3d,Knots);
1236                 BSplCLib::Reparametrize(fC0,lC0,Knots);
1237                 bs2d->SetKnots(Knots);
1238                 GAC2d.Load(bs2d,f3d,l3d);
1239                 curPC = bs2d;
1240                 Standard_Boolean updatepcsov = updatepc;
1241                 updatepc = Standard_True;
1242
1243                 Standard_Real error1 = ComputeTol(HC, HC2d, HS, NCONTROL);
1244                 if(error1 > error) {
1245                   bs2d = bs2dsov;
1246                   GAC2d.Load(bs2d,f3d,l3d);
1247                   curPC = bs2d;
1248                   updatepc = updatepcsov;
1249                   isANA = Standard_True;
1250                 }
1251                 else {
1252                   error = error1;
1253                 }
1254               }
1255
1256               //check, if new BSpline "good" or not --------- IFV, Jan of 2000
1257               GeomAbs_Shape cont = bs2d->Continuity();
1258               Standard_Boolean IsBad = Standard_False;
1259
1260               if(cont > GeomAbs_C0 && error > Max(1.e-3,Tolerance)) {
1261                 Standard_Integer NbKnots = bs2d->NbKnots();
1262                 TColStd_Array1OfReal Knots(1,NbKnots);
1263                 bs2d->Knots(Knots);
1264                 Standard_Real critratio = 10.; 
1265                 Standard_Real dtprev = Knots(2) - Knots(1), dtratio = 1.;
1266                 Standard_Real dtmin = dtprev;
1267                 Standard_Real dtcur;
1268                 for(Standard_Integer j = 2; j < NbKnots; j++) {
1269                   dtcur = Knots(j+1) - Knots(j);
1270                   dtmin = Min(dtmin, dtcur);
1271
1272                   if(IsBad) continue;
1273
1274                   if(dtcur > dtprev) dtratio = dtcur/dtprev;
1275                   else dtratio = dtprev/dtcur;
1276                   if(dtratio > critratio) {IsBad = Standard_True;}
1277                   dtprev = dtcur;
1278
1279                 }
1280                 if(IsBad) {
1281                   // To avoid failures in Approx_CurvilinearParameter 
1282                   bs2d->Resolution(Max(1.e-3,Tolerance), dtcur);
1283                   if(dtmin < dtcur) IsBad = Standard_False;
1284                 }
1285               }
1286
1287
1288               if(IsBad ) { //if BSpline "bad", try to reparametrize it
1289                 // by its curve length
1290
1291                 //            GeomAbs_Shape cont = bs2d->Continuity();
1292                 if(cont > GeomAbs_C2) cont = GeomAbs_C2;
1293                 Standard_Integer maxdeg = bs2d->Degree();
1294                 if(maxdeg == 1) maxdeg = 14;
1295                 Approx_CurvilinearParameter AppCurPar(HC2d, HS, Max(1.e-3,Tolerance),
1296                   cont, maxdeg, 10);
1297                 if(AppCurPar.IsDone() || AppCurPar.HasResult()) {
1298                   bs2d = AppCurPar.Curve2d1();
1299                   GAC2d.Load(bs2d,f3d,l3d);
1300                   curPC = bs2d;
1301
1302                   if(Abs(bs2d->FirstParameter() - fC0) > Tol2d ||
1303                     Abs(bs2d->LastParameter() - lC0) > Tol2d    ) {
1304                       Standard_Integer NbKnots = bs2d->NbKnots();
1305                       TColStd_Array1OfReal Knots(1,NbKnots);
1306                       bs2d->Knots(Knots);
1307                       //                  BSplCLib::Reparametrize(f3d,l3d,Knots);
1308                       BSplCLib::Reparametrize(fC0,lC0,Knots);
1309                       bs2d->SetKnots(Knots);
1310                       GAC2d.Load(bs2d,f3d,l3d);
1311                       curPC = bs2d;
1312
1313                   }
1314                 }
1315               }
1316
1317
1318             }
1319         }
1320
1321
1322         if(goodpc){
1323           //      Approx_SameParameter SameP(HC,HC2d,HS,Tolerance);
1324           Standard_Real aTol = (isANA && isBSP) ? 1.e-7 : Tolerance;
1325           const Handle(Adaptor3d_HCurve)& aHCurv = HC; // to avoid ambiguity
1326           const Handle(Adaptor2d_HCurve2d)& aHCurv2d = HC2d; // to avoid ambiguity
1327           Approx_SameParameter SameP(aHCurv,aHCurv2d,HS,aTol);
1328
1329           if (SameP.IsSameParameter()) {
1330             maxdist = Max(maxdist,SameP.TolReached());
1331             if(updatepc){
1332               if (i == 0) GCurve->PCurve(curPC);
1333               else GCurve->PCurve2(curPC);
1334             }
1335           }
1336           else if (SameP.IsDone()) {
1337             Standard_Real tolreached = SameP.TolReached();
1338             if(tolreached <= error) {
1339               curPC = SameP.Curve2d();
1340               updatepc = Standard_True;
1341               maxdist = Max(maxdist,tolreached);
1342             }
1343             else {
1344               maxdist = Max(maxdist, error);
1345             }
1346             if(updatepc){
1347               if (i == 0) GCurve->PCurve(curPC);
1348               else GCurve->PCurve2(curPC);
1349             }
1350           }
1351           else
1352           {
1353             //Approx_SameParameter has failed.
1354             //Consequently, the situation might be,
1355             //when 3D and 2D-curve do not have same-range.
1356             GeomLib::SameRange( Tol2d, PC[i], 
1357                                 GCurve->First(), GCurve->Last(),
1358                                 f3d,l3d,curPC);
1359
1360             if (i == 0) GCurve->PCurve(curPC);
1361             else GCurve->PCurve2(curPC);
1362
1363             IsSameP = 0;
1364           }
1365
1366         }
1367         else IsSameP = 0;
1368
1369         //  Modified by skv - Thu Jun  3 12:39:19 2004 OCC5898 Begin
1370         if (!IsSameP) {
1371           if (anEdgeTol >= error) {
1372             maxdist = Max(maxdist, anEdgeTol);
1373             IsSameP = Standard_True;
1374           }
1375         }
1376         //  Modified by skv - Thu Jun  3 12:39:20 2004 OCC5898 End
1377       }
1378     }
1379     It.Next() ;
1380   }
1381   B.Range(AnEdge,f3d,l3d);
1382   B.SameRange(AnEdge,Standard_True);
1383   if ( IsSameP) {
1384     // Reduce eventually the tolerance of the edge, as
1385     // all its representations are processed (except for some associated
1386     // to planes and not stored in the edge !) 
1387     // The same cannot be done with vertices that cannot be enlarged 
1388     // or left as is.
1389     if (YaPCu) {
1390       // Avoid setting too small tolerances.
1391       maxdist = Max(maxdist,Precision::Confusion());
1392       TopoDS_Vertex V1,V2;
1393       TopExp::Vertices(AnEdge,V1,V2);
1394       if (!V1.IsNull())
1395         B.UpdateVertex(V1,maxdist);
1396       if (!V2.IsNull())
1397         B.UpdateVertex(V2,maxdist);
1398       TE->Modified(Standard_True);
1399       TE->Tolerance(maxdist);
1400     }
1401     B.SameParameter(AnEdge,Standard_True);
1402   }
1403 }
1404
1405 //=======================================================================
1406 //function : UpdateTolerances
1407 //purpose  : 
1408 //=======================================================================
1409 void  BRepLib::UpdateTolerances(const TopoDS_Shape& aShape,
1410   const Standard_Boolean verifyTolerance) 
1411 {
1412
1413   // Harmonize tolerances
1414   // with rule Tolerance(VERTEX)>=Tolerance(EDGE)>=Tolerance(FACE)
1415   BRep_Builder B;
1416   Standard_Real tol=0;
1417   if (verifyTolerance) {
1418     // Set tolerance to its minimum value
1419     Handle(Geom_Surface) S;
1420     TopLoc_Location l;
1421     TopExp_Explorer ex;
1422     Bnd_Box aB;
1423     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, dMax;
1424     for (ex.Init(aShape, TopAbs_FACE); ex.More(); ex.Next()) {
1425       const TopoDS_Face& curf=TopoDS::Face(ex.Current());
1426       S = BRep_Tool::Surface(curf, l);
1427       if (!S.IsNull()) {
1428         aB.SetVoid();
1429         BRepBndLib::Add(curf,aB);
1430         if (S->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
1431           S = Handle(Geom_RectangularTrimmedSurface)::DownCast (S)->BasisSurface();
1432         }
1433         GeomAdaptor_Surface AS(S);
1434         switch (AS.GetType()) {
1435         case GeomAbs_Plane: 
1436         case GeomAbs_Cylinder: 
1437         case GeomAbs_Cone: 
1438           {
1439             tol=Precision::Confusion();
1440             break;
1441           }
1442         case GeomAbs_Sphere: 
1443         case GeomAbs_Torus: 
1444           {
1445             tol=Precision::Confusion()*2;
1446             break;
1447           }
1448         default:
1449           tol=Precision::Confusion()*4;
1450         }
1451         if (!aB.IsWhole()) {
1452           aB.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
1453           dMax=1.;
1454           if (!aB.IsOpenXmin() && !aB.IsOpenXmax()) dMax=aXmax-aXmin;
1455           if (!aB.IsOpenYmin() && !aB.IsOpenYmax()) aYmin=aYmax-aYmin;
1456           if (!aB.IsOpenZmin() && !aB.IsOpenZmax()) aZmin=aZmax-aZmin;
1457           if (aYmin>dMax) dMax=aYmin;
1458           if (aZmin>dMax) dMax=aZmin;
1459           tol=tol*dMax;
1460           // Do not process tolerances > 1.
1461           if (tol>1.) tol=0.99;
1462         }
1463         const Handle(BRep_TFace)& Tf = *((Handle(BRep_TFace)*)&curf.TShape());
1464         Tf->Tolerance(tol);
1465       }
1466     }
1467   }
1468
1469   //Process edges
1470   TopTools_IndexedDataMapOfShapeListOfShape parents;
1471   TopExp::MapShapesAndAncestors(aShape, TopAbs_EDGE, TopAbs_FACE, parents);
1472   TopTools_ListIteratorOfListOfShape lConx;
1473   Standard_Integer iCur;
1474   for (iCur=1; iCur<=parents.Extent(); iCur++) {
1475     tol=0;
1476     for (lConx.Initialize(parents(iCur)); lConx.More(); lConx.Next()) {
1477       tol=Max(tol, BRep_Tool::Tolerance(TopoDS::Face(lConx.Value())));
1478     }
1479     // Update can only increase tolerance, so if the edge has a greater
1480     //  tolerance than its faces it is not concerned
1481     B.UpdateEdge(TopoDS::Edge(parents.FindKey(iCur)), tol);
1482   }
1483
1484   //Vertices are processed
1485   const Standard_Real BigTol = 1.e10;
1486   parents.Clear();
1487   TopExp::MapShapesAndUniqueAncestors(aShape, TopAbs_VERTEX, TopAbs_EDGE, parents);
1488   TColStd_MapOfTransient Initialized;
1489   Standard_Integer nbV = parents.Extent();
1490   for (iCur=1; iCur<=nbV; iCur++) {
1491     tol=0;
1492     const TopoDS_Vertex& V = TopoDS::Vertex(parents.FindKey(iCur));
1493     Bnd_Box box;
1494     box.Add(BRep_Tool::Pnt(V));
1495     gp_Pnt p3d;
1496     for (lConx.Initialize(parents(iCur)); lConx.More(); lConx.Next()) {
1497       const TopoDS_Edge& E = TopoDS::Edge(lConx.Value());
1498       tol=Max(tol, BRep_Tool::Tolerance(E));
1499       if(tol > BigTol) continue;
1500       if(!BRep_Tool::SameRange(E)) continue;
1501       Standard_Real par = BRep_Tool::Parameter(V,E);
1502       Handle(BRep_TEdge)& TE = *((Handle(BRep_TEdge)*)&E.TShape());
1503       BRep_ListIteratorOfListOfCurveRepresentation itcr(TE->Curves());
1504       const TopLoc_Location& Eloc = E.Location();
1505       while (itcr.More()) {
1506         // For each CurveRepresentation, check the provided parameter
1507         const Handle(BRep_CurveRepresentation)& cr = itcr.Value();
1508         const TopLoc_Location& loc = cr->Location();
1509         TopLoc_Location L = (Eloc * loc);
1510         if (cr->IsCurve3D()) {
1511           const Handle(Geom_Curve)& C = cr->Curve3D();
1512           if (!C.IsNull()) { // edge non degenerated
1513             p3d = C->Value(par);
1514             p3d.Transform(L.Transformation());
1515             box.Add(p3d);
1516           }
1517         }
1518         else if (cr->IsCurveOnSurface()) {
1519           const Handle(Geom_Surface)& Su = cr->Surface();
1520           const Handle(Geom2d_Curve)& PC = cr->PCurve();
1521           Handle(Geom2d_Curve) PC2;
1522           if (cr->IsCurveOnClosedSurface()) {
1523             PC2 = cr->PCurve2();
1524           }
1525           gp_Pnt2d p2d = PC->Value(par);
1526           p3d = Su->Value(p2d.X(),p2d.Y());
1527           p3d.Transform(L.Transformation());
1528           box.Add(p3d);
1529           if (!PC2.IsNull()) {
1530             p2d = PC2->Value(par);
1531             p3d = Su->Value(p2d.X(),p2d.Y());
1532             p3d.Transform(L.Transformation());
1533             box.Add(p3d);
1534           }
1535         }
1536         itcr.Next();
1537       }
1538     }
1539     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
1540     box.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
1541     aXmax -= aXmin; aYmax -= aYmin; aZmax -= aZmin;
1542     tol = Max(tol,sqrt(aXmax*aXmax+aYmax*aYmax+aZmax*aZmax));
1543     tol += 2.*Epsilon(tol);
1544     if (verifyTolerance) {
1545       // ASet minimum value of the tolerance 
1546       // Attention to sharing of the vertex by other shapes
1547       const Handle(BRep_TVertex)& TV = *((Handle(BRep_TVertex)*)&V.TShape());
1548       if (Initialized.Add(TV)) 
1549         TV->Tolerance(tol);
1550       else 
1551         B.UpdateVertex(V, tol);
1552     }
1553     else {
1554       // Update can only increase tolerance, so if the edge has a greater
1555       //  tolerance than its faces it is not concerned
1556       B.UpdateVertex(V, tol);
1557     }
1558   }
1559 }
1560
1561 //=======================================================================
1562 //function : UpdateInnerTolerances
1563 //purpose  : 
1564 //=======================================================================
1565 void  BRepLib::UpdateInnerTolerances(const TopoDS_Shape& aShape)
1566 {
1567   TopTools_IndexedDataMapOfShapeListOfShape EFmap;
1568   TopExp::MapShapesAndAncestors(aShape, TopAbs_EDGE, TopAbs_FACE, EFmap);
1569   BRep_Builder BB;
1570   for (Standard_Integer i = 1; i <= EFmap.Extent(); i++)
1571   {
1572     TopoDS_Edge anEdge = TopoDS::Edge(EFmap.FindKey(i));
1573     TopoDS_Vertex V1, V2;
1574     TopExp::Vertices(anEdge, V1, V2);
1575     Standard_Real fpar, lpar;
1576     BRep_Tool::Range(anEdge, fpar, lpar);
1577     Standard_Real TolEdge = BRep_Tool::Tolerance(anEdge);
1578     gp_Pnt Pnt1, Pnt2;
1579     Handle(BRepAdaptor_HCurve) anHCurve = new BRepAdaptor_HCurve();
1580     anHCurve->ChangeCurve().Initialize(anEdge);
1581     if (!V1.IsNull())
1582       Pnt1 = BRep_Tool::Pnt(V1);
1583     if (!V2.IsNull())
1584       Pnt2 = BRep_Tool::Pnt(V2);
1585     
1586     if (!BRep_Tool::Degenerated(anEdge) &&
1587         EFmap(i).Extent() > 0)
1588     {
1589       NCollection_Sequence<Handle(Adaptor3d_HCurve)> theRep;
1590       theRep.Append(anHCurve);
1591       TopTools_ListIteratorOfListOfShape itl(EFmap(i));
1592       for (; itl.More(); itl.Next())
1593       {
1594         const TopoDS_Face& aFace = TopoDS::Face(itl.Value());
1595         Handle(BRepAdaptor_HCurve) anHCurvOnSurf = new BRepAdaptor_HCurve();
1596         anHCurvOnSurf->ChangeCurve().Initialize(anEdge, aFace);
1597         theRep.Append(anHCurvOnSurf);
1598       }
1599       
1600       const Standard_Integer NbSamples = (BRep_Tool::SameParameter(anEdge))? 23 : 2;
1601       Standard_Real delta = (lpar - fpar)/(NbSamples-1);
1602       Standard_Real MaxDist = 0.;
1603       for (Standard_Integer j = 2; j <= theRep.Length(); j++)
1604       {
1605         for (Standard_Integer k = 0; k <= NbSamples; k++)
1606         {
1607           Standard_Real ParamOnCenter = (k == NbSamples)? lpar :
1608             fpar + k*delta;
1609           gp_Pnt Center = theRep(1)->Value(ParamOnCenter);
1610           Standard_Real ParamOnCurve = (BRep_Tool::SameParameter(anEdge))? ParamOnCenter
1611             : ((k == 0)? theRep(j)->FirstParameter() : theRep(j)->LastParameter());
1612           gp_Pnt aPoint = theRep(j)->Value(ParamOnCurve);
1613           Standard_Real aDist = Center.Distance(aPoint);
1614           //aDist *= 1.1;
1615           aDist += 2.*Epsilon(aDist);
1616           if (aDist > MaxDist)
1617             MaxDist = aDist;
1618
1619           //Update tolerances of vertices
1620           if (k == 0 && !V1.IsNull())
1621           {
1622             Standard_Real aDist1 = Pnt1.Distance(aPoint);
1623             aDist1 += 2.*Epsilon(aDist1);
1624             BB.UpdateVertex(V1, aDist1);
1625           }
1626           if (k == NbSamples && !V2.IsNull())
1627           {
1628             Standard_Real aDist2 = Pnt2.Distance(aPoint);
1629             aDist2 += 2.*Epsilon(aDist2);
1630             BB.UpdateVertex(V2, aDist2);
1631           }
1632         }
1633       }
1634       BB.UpdateEdge(anEdge, MaxDist);
1635     }
1636     TolEdge = BRep_Tool::Tolerance(anEdge);
1637     if (!V1.IsNull())
1638     {
1639       gp_Pnt End1 = anHCurve->Value(fpar);
1640       Standard_Real dist1 = Pnt1.Distance(End1);
1641       dist1 += 2.*Epsilon(dist1);
1642       BB.UpdateVertex(V1, Max(dist1, TolEdge));
1643     }
1644     if (!V2.IsNull())
1645     {
1646       gp_Pnt End2 = anHCurve->Value(lpar);
1647       Standard_Real dist2 = Pnt2.Distance(End2);
1648       dist2 += 2.*Epsilon(dist2);
1649       BB.UpdateVertex(V2, Max(dist2, TolEdge));
1650     }
1651   }
1652 }
1653
1654 //=======================================================================
1655 //function : OrientClosedSolid
1656 //purpose  : 
1657 //=======================================================================
1658 Standard_Boolean BRepLib::OrientClosedSolid(TopoDS_Solid& solid) 
1659 {
1660   // Set material inside the solid
1661   BRepClass3d_SolidClassifier where(solid);
1662   where.PerformInfinitePoint(Precision::Confusion());
1663   if (where.State()==TopAbs_IN) {
1664     solid.Reverse();
1665   }
1666   else if (where.State()==TopAbs_ON || where.State()==TopAbs_UNKNOWN) 
1667     return Standard_False;
1668
1669   return Standard_True;
1670 }
1671
1672 // Structure for calculation of properties, necessary for decision about continuity
1673 class SurfaceProperties
1674 {
1675 public:
1676   SurfaceProperties(const Handle(Geom_Surface)& theSurface,
1677                     const gp_Trsf&              theSurfaceTrsf,
1678                     const Handle(Geom2d_Curve)& theCurve2D,
1679                     const Standard_Boolean      theReversed)
1680     : mySurfaceProps(theSurface, 2, Precision::Confusion()),
1681       mySurfaceTrsf(theSurfaceTrsf),
1682       myCurve2d(theCurve2D),
1683       myIsReversed(theReversed)
1684   {}
1685
1686   // Calculate derivatives on surface related to the point on curve
1687   void Calculate(const Standard_Real theParamOnCurve)
1688   {
1689     gp_Pnt2d aUV;
1690     myCurve2d->D1(theParamOnCurve, aUV, myCurveTangent);
1691     mySurfaceProps.SetParameters(aUV.X(), aUV.Y());
1692   }
1693
1694   // Returns point just calculated
1695   gp_Pnt Value() 
1696   { return mySurfaceProps.Value().Transformed(mySurfaceTrsf); }
1697
1698   // Calculate a derivative orthogonal to curve's tangent vector
1699   gp_Vec Derivative()
1700   {
1701     gp_Vec aDeriv;
1702     // direction orthogonal to tangent vector of the curve
1703     gp_Vec2d anOrtho(-myCurveTangent.Y(), myCurveTangent.X());
1704     Standard_Real aLen = anOrtho.Magnitude();
1705     if (aLen < Precision::Confusion())
1706       return aDeriv;
1707     anOrtho /= aLen;
1708     if (myIsReversed)
1709       anOrtho.Reverse();
1710
1711     aDeriv.SetLinearForm(anOrtho.X(), mySurfaceProps.D1U(),
1712                          anOrtho.Y(), mySurfaceProps.D1V());
1713     return aDeriv.Transformed(mySurfaceTrsf);
1714   }
1715
1716   // Calculate principal curvatures, which consist of minimal and maximal normal curvatures and
1717   // the directions on the tangent plane (principal direction) where the extremums are reached
1718   void Curvature(gp_Dir& thePrincipalDir1, Standard_Real& theCurvature1,
1719                  gp_Dir& thePrincipalDir2, Standard_Real& theCurvature2)
1720   {
1721     mySurfaceProps.CurvatureDirections(thePrincipalDir1, thePrincipalDir2);
1722     theCurvature1 = mySurfaceProps.MaxCurvature();
1723     theCurvature2 = mySurfaceProps.MinCurvature();
1724     if (myIsReversed)
1725     {
1726       theCurvature1 = -theCurvature1;
1727       theCurvature2 = -theCurvature2;
1728     }
1729     thePrincipalDir1.Transform(mySurfaceTrsf);
1730     thePrincipalDir2.Transform(mySurfaceTrsf);
1731   }
1732
1733 private:
1734   GeomLProp_SLProps    mySurfaceProps; // properties calculator
1735   gp_Trsf              mySurfaceTrsf;
1736   Handle(Geom2d_Curve) myCurve2d;
1737   Standard_Boolean     myIsReversed; // the face based on the surface is reversed
1738
1739   // tangent vector to Pcurve in UV
1740   gp_Vec2d myCurveTangent;
1741 };
1742
1743 //=======================================================================
1744 //function : tgtfaces
1745 //purpose  : check the angle at the border between two squares.
1746 //           Two shares should have a shared front edge.
1747 //=======================================================================
1748 static GeomAbs_Shape tgtfaces(const TopoDS_Edge& Ed,
1749                               const TopoDS_Face& F1,
1750                               const TopoDS_Face& F2,
1751                               const Standard_Real theAngleTol)
1752 {
1753   Standard_Boolean isSeam = F1.IsEqual(F2);
1754
1755   TopoDS_Edge E = Ed;
1756
1757   // Check if pcurves exist on both faces of edge
1758   Standard_Real aFirst,aLast;
1759   E.Orientation(TopAbs_FORWARD);
1760   Handle(Geom2d_Curve) aCurve1 = BRep_Tool::CurveOnSurface(E, F1, aFirst, aLast);
1761   if(aCurve1.IsNull())
1762     return GeomAbs_C0;
1763
1764   if (isSeam)
1765     E.Orientation(TopAbs_REVERSED);
1766   Handle(Geom2d_Curve) aCurve2 = BRep_Tool::CurveOnSurface(E, F2, aFirst, aLast);
1767   if(aCurve2.IsNull())
1768     return GeomAbs_C0;
1769
1770   TopLoc_Location aLoc1, aLoc2;
1771   Handle(Geom_Surface) aSurface1 = BRep_Tool::Surface(F1, aLoc1);
1772   const gp_Trsf& aSurf1Trsf = aLoc1.Transformation();
1773   Handle(Geom_Surface) aSurface2 = BRep_Tool::Surface(F2, aLoc2);
1774   const gp_Trsf& aSurf2Trsf = aLoc2.Transformation();
1775
1776   if (aSurface1->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
1777     aSurface1 = Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface1)->BasisSurface();
1778   if (aSurface2->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
1779     aSurface2 = Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface2)->BasisSurface();
1780
1781   // seam edge on elementary surface is always CN
1782   Standard_Boolean isElementary =
1783     (aSurface1->IsKind(STANDARD_TYPE(Geom_ElementarySurface)) &&
1784      aSurface2->IsKind(STANDARD_TYPE(Geom_ElementarySurface)));
1785   if (isSeam && isElementary)
1786   {
1787     return GeomAbs_CN;
1788   }
1789
1790   SurfaceProperties aSP1(aSurface1, aSurf1Trsf, aCurve1, F1.Orientation() == TopAbs_REVERSED);
1791   SurfaceProperties aSP2(aSurface2, aSurf2Trsf, aCurve2, F2.Orientation() == TopAbs_REVERSED);
1792
1793   Standard_Real f, l, eps;
1794   BRep_Tool::Range(E,f,l);
1795   Extrema_LocateExtPC ext;
1796   Handle(BRepAdaptor_HCurve) aHC2;
1797
1798   eps = (l - f)/100.;
1799   f += eps; // to avoid calculations on  
1800   l -= eps; // points of pointed squares.
1801
1802   const Standard_Real anAngleTol2 = theAngleTol * theAngleTol;
1803
1804   gp_Vec aDer1, aDer2;
1805   gp_Vec aNorm1;
1806   Standard_Real aSqLen1, aSqLen2;
1807   gp_Dir aCrvDir1[2], aCrvDir2[2];
1808   Standard_Real aCrvLen1[2], aCrvLen2[2];
1809
1810   GeomAbs_Shape aCont = (isElementary ? GeomAbs_CN : GeomAbs_C2);
1811   GeomAbs_Shape aCurCont;
1812   Standard_Real u;
1813   for (Standard_Integer i = 0; i <= 20 && aCont > GeomAbs_C0; i++)
1814   {
1815     // First suppose that this is sameParameter
1816     u = f + (l-f)*i/20;
1817
1818     // Check conditions for G1 and C1 continuity:
1819     // * calculate a derivative in tangent plane of each surface
1820     //   orthogonal to curve's tangent vector
1821     // * continuity is C1 if the vectors are equal
1822     // * continuity is G1 if the vectors are just parallel
1823     aCurCont = GeomAbs_C0;
1824
1825     aSP1.Calculate(u);
1826     aSP2.Calculate(u);
1827
1828     aDer1 = aSP1.Derivative();
1829     aSqLen1 = aDer1.SquareMagnitude();
1830     aDer2 = aSP2.Derivative();
1831     aSqLen2 = aDer2.SquareMagnitude();
1832     Standard_Boolean isSmoothSuspect = (aDer1.CrossSquareMagnitude(aDer2) <= anAngleTol2 * aSqLen1 * aSqLen2);
1833     if (!isSmoothSuspect)
1834     {
1835       // Refine by projection
1836       if (aHC2.IsNull())
1837       {
1838         // adaptor for pcurve on the second surface
1839         aHC2 = new BRepAdaptor_HCurve(BRepAdaptor_Curve(E, F2));
1840         ext.Initialize(aHC2->Curve(), f, l, Precision::PConfusion());
1841       }
1842       ext.Perform(aSP1.Value(), u);
1843       if (ext.IsDone() && ext.IsMin())
1844       {
1845         const Extrema_POnCurv& poc = ext.Point();
1846         aSP2.Calculate(poc.Parameter());
1847         aDer2 = aSP2.Derivative();
1848         aSqLen2 = aDer2.SquareMagnitude();
1849       }
1850       isSmoothSuspect = (aDer1.CrossSquareMagnitude(aDer2) <= anAngleTol2 * aSqLen1 * aSqLen2);
1851     }
1852     if (isSmoothSuspect)
1853     {
1854       aCurCont = GeomAbs_G1;
1855       if (Abs(Sqrt(aSqLen1) - Sqrt(aSqLen2)) < Precision::Confusion() &&
1856           aDer1.Dot(aDer2) > Precision::SquareConfusion()) // <= check vectors are codirectional
1857         aCurCont = GeomAbs_C1;
1858     }
1859     else
1860       return GeomAbs_C0;
1861
1862     if (aCont < GeomAbs_G2)
1863       continue; // no need further processing, because maximal continuity is less than G2
1864
1865     // Check conditions for G2 and C2 continuity:
1866     // * calculate principal curvatures on each surface
1867     // * continuity is C2 if directions of principal curvatures are equal on differenct surfaces
1868     // * continuity is G2 if directions of principal curvatures are just parallel
1869     //   and values of curvatures are the same
1870     aSP1.Curvature(aCrvDir1[0], aCrvLen1[0], aCrvDir1[1], aCrvLen1[1]);
1871     aSP2.Curvature(aCrvDir2[0], aCrvLen2[0], aCrvDir2[1], aCrvLen2[1]);
1872     for (Standard_Integer aStep = 0; aStep <= 1; ++aStep)
1873     {
1874       if (aCrvDir1[0].XYZ().CrossSquareMagnitude(aCrvDir2[aStep].XYZ()) <= Precision::SquareConfusion() &&
1875           Abs(aCrvLen1[0] - aCrvLen2[aStep]) < Precision::Confusion() &&
1876           aCrvDir1[1].XYZ().CrossSquareMagnitude(aCrvDir2[1 - aStep].XYZ()) <= Precision::SquareConfusion() &&
1877           Abs(aCrvLen1[1] - aCrvLen2[1 - aStep]) < Precision::Confusion())
1878       {
1879         if (aCurCont == GeomAbs_C1 &&
1880             aCrvDir1[0].Dot(aCrvDir2[aStep]) > Precision::Confusion() &&
1881             aCrvDir1[1].Dot(aCrvDir2[1 - aStep]) > Precision::Confusion())
1882           aCurCont = GeomAbs_C2;
1883         else
1884           aCurCont = GeomAbs_G2;
1885         break;
1886       }
1887     }
1888
1889     if (aCurCont < aCont)
1890       aCont = aCurCont;
1891   }
1892
1893   // according to the list of supported elementary surfaces,
1894   // if the continuity is C2, than it is totally CN
1895   if (isElementary && aCont == GeomAbs_C2)
1896     aCont = GeomAbs_CN;
1897   return aCont;
1898 }
1899
1900 //=======================================================================
1901 // function : EncodeRegularity
1902 // purpose  : Code the regularities on all edges of the shape, boundary of 
1903 //            two faces that do not have it.
1904 //            Takes into account that compound may consists of same solid
1905 //            placed with different transformations
1906 //=======================================================================
1907 static void EncodeRegularity(const TopoDS_Shape&        theShape,
1908                              const Standard_Real        theTolAng,
1909                              TopTools_MapOfShape&       theMap,
1910                              const TopTools_MapOfShape& theEdgesToEncode = TopTools_MapOfShape())
1911 {
1912   TopoDS_Shape aShape = theShape;
1913   TopLoc_Location aNullLoc;
1914   aShape.Location(aNullLoc); // nullify location
1915   if (!theMap.Add(aShape))
1916     return; // do not need to process shape twice
1917
1918   if (aShape.ShapeType() == TopAbs_COMPOUND ||
1919       aShape.ShapeType() == TopAbs_COMPSOLID)
1920   {
1921     for (TopoDS_Iterator it(aShape); it.More(); it.Next())
1922       EncodeRegularity(it.Value(), theTolAng, theMap, theEdgesToEncode);
1923     return;
1924   }
1925
1926   try {
1927     OCC_CATCH_SIGNALS
1928
1929     TopTools_IndexedDataMapOfShapeListOfShape M;
1930     TopExp::MapShapesAndAncestors(aShape, TopAbs_EDGE, TopAbs_FACE, M);
1931     TopTools_ListIteratorOfListOfShape It;
1932     TopExp_Explorer Ex;
1933     TopoDS_Face F1,F2;
1934     Standard_Boolean found;
1935     for (Standard_Integer i = 1; i <= M.Extent(); i++){
1936       TopoDS_Edge E = TopoDS::Edge(M.FindKey(i));
1937       if (!theEdgesToEncode.IsEmpty())
1938       {
1939         // process only the edges from the list to update their regularity
1940         TopoDS_Shape aPureEdge = E.Located(aNullLoc);
1941         aPureEdge.Orientation(TopAbs_FORWARD);
1942         if (!theEdgesToEncode.Contains(aPureEdge))
1943           continue;
1944       }
1945
1946       found = Standard_False;                                     
1947       F1.Nullify();
1948       for (It.Initialize(M.FindFromIndex(i)); It.More() && !found; It.Next()){
1949         if (F1.IsNull()) { F1 = TopoDS::Face(It.Value()); }
1950         else {
1951           const TopoDS_Face& aTmpF2 = TopoDS::Face(It.Value());
1952           if (!F1.IsSame(aTmpF2)){
1953             found = Standard_True;
1954             F2 = aTmpF2;
1955           }
1956         }
1957       }
1958       if (!found && !F1.IsNull()){//is it a sewing edge?
1959         TopAbs_Orientation orE = E.Orientation();
1960         TopoDS_Edge curE;
1961         for (Ex.Init(F1, TopAbs_EDGE); Ex.More() && !found; Ex.Next()){
1962           curE = TopoDS::Edge(Ex.Current());
1963           if (E.IsSame(curE) && orE != curE.Orientation()) {
1964             found = Standard_True;
1965             F2 = F1;
1966           }
1967         }
1968       }
1969       if (found)
1970         BRepLib::EncodeRegularity(E, F1, F2, theTolAng);
1971     }
1972   }
1973   catch (Standard_Failure const& anException) {
1974 #ifdef OCCT_DEBUG
1975     cout << "Warning: Exception in BRepLib::EncodeRegularity(): ";
1976     anException.Print(cout);
1977     cout << endl;
1978 #endif
1979     (void)anException;
1980   }
1981 }
1982
1983 //=======================================================================
1984 // function : EncodeRegularity
1985 // purpose  : code the regularities on all edges of the shape, boundary of 
1986 //            two faces that do not have it.
1987 //=======================================================================
1988
1989 void BRepLib::EncodeRegularity(const TopoDS_Shape& S,
1990   const Standard_Real TolAng)
1991 {
1992   TopTools_MapOfShape aMap;
1993   ::EncodeRegularity(S, TolAng, aMap);
1994 }
1995
1996 //=======================================================================
1997 // function : EncodeRegularity
1998 // purpose  : code the regularities on all edges in the list that do not 
1999 //            have it, and which are boundary of two faces on the shape.
2000 //=======================================================================
2001
2002 void BRepLib::EncodeRegularity(const TopoDS_Shape& S,
2003   const TopTools_ListOfShape& LE,
2004   const Standard_Real TolAng)
2005 {
2006   // Collect edges without location and orientation
2007   TopTools_MapOfShape aPureEdges;
2008   TopLoc_Location aNullLoc;
2009   TopTools_ListIteratorOfListOfShape anEdgeIt(LE);
2010   for (; anEdgeIt.More(); anEdgeIt.Next())
2011   {
2012     TopoDS_Shape anEdge = anEdgeIt.Value();
2013     anEdge.Location(aNullLoc);
2014     anEdge.Orientation(TopAbs_FORWARD);
2015     aPureEdges.Add(anEdge);
2016   }
2017
2018   TopTools_MapOfShape aMap;
2019   ::EncodeRegularity(S, TolAng, aMap, aPureEdges);
2020 }
2021
2022 //=======================================================================
2023 // function : EncodeRegularity
2024 // purpose  : code the regularity between 2 faces connected by edge 
2025 //=======================================================================
2026
2027 void BRepLib::EncodeRegularity(TopoDS_Edge& E,
2028   const TopoDS_Face& F1,
2029   const TopoDS_Face& F2,
2030   const Standard_Real TolAng)
2031 {
2032   BRep_Builder B;
2033   if(BRep_Tool::Continuity(E,F1,F2)<=GeomAbs_C0){
2034     try {
2035       GeomAbs_Shape aCont = tgtfaces(E, F1, F2, TolAng);
2036       B.Continuity(E,F1,F2,aCont);
2037       
2038     }
2039     catch(Standard_Failure)
2040     {
2041 #ifdef OCCT_DEBUG
2042       cout << "Failure: Exception in BRepLib::EncodeRegularity" << endl;
2043 #endif
2044     }
2045   }
2046 }
2047
2048 //=======================================================================
2049 // function : EnsureNormalConsistency
2050 // purpose  : Corrects the normals in Poly_Triangulation of faces.
2051 //              Returns TRUE if any correction is done.
2052 //=======================================================================
2053 Standard_Boolean BRepLib::
2054             EnsureNormalConsistency(const TopoDS_Shape& theShape,
2055                                     const Standard_Real theAngTol,
2056                                     const Standard_Boolean theForceComputeNormals)
2057 {
2058   const Standard_Real aThresDot = cos(theAngTol);
2059
2060   Standard_Boolean aRetVal = Standard_False, isNormalsFound = Standard_False;
2061
2062   // compute normals if they are absent
2063   TopExp_Explorer anExpFace(theShape,TopAbs_FACE);
2064   for (; anExpFace.More(); anExpFace.Next())
2065   {
2066     const TopoDS_Face& aFace = TopoDS::Face(anExpFace.Current());
2067     const Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
2068     if(aSurf.IsNull())
2069       continue;
2070     TopLoc_Location aLoc;
2071     const Handle(Poly_Triangulation)& aPT = BRep_Tool::Triangulation(aFace, aLoc);
2072     if(aPT.IsNull())
2073       continue;
2074     if (!theForceComputeNormals && aPT->HasNormals())
2075     {
2076       isNormalsFound = Standard_True;
2077       continue;
2078     }
2079
2080     GeomLProp_SLProps aSLP(aSurf, 2, Precision::Confusion());
2081     const Standard_Integer anArrDim = 3*aPT->NbNodes();
2082     Handle(TShort_HArray1OfShortReal) aNormArr = new TShort_HArray1OfShortReal(1, anArrDim);
2083     Standard_Integer anNormInd = aNormArr->Lower();
2084     for(Standard_Integer i = aPT->UVNodes().Lower(); i <= aPT->UVNodes().Upper(); i++)
2085     {
2086       const gp_Pnt2d &aP2d = aPT->UVNodes().Value(i);
2087       aSLP.SetParameters(aP2d.X(), aP2d.Y());
2088
2089       gp_XYZ aNorm(0.,0.,0.);
2090       if(!aSLP.IsNormalDefined())
2091       {
2092 #ifdef OCCT_DEBUG
2093         cout << "BRepLib::EnsureNormalConsistency(): Cannot find normal!" << endl;
2094 #endif
2095       }
2096       else
2097       {
2098         aNorm = aSLP.Normal().XYZ();
2099         if (aFace.Orientation() == TopAbs_REVERSED)
2100           aNorm.Reverse();
2101       }
2102       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.X());
2103       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.Y());
2104       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.Z());
2105     }
2106
2107     aRetVal = Standard_True;
2108     isNormalsFound = Standard_True;
2109     aPT->SetNormals(aNormArr);
2110   }
2111
2112   if(!isNormalsFound)
2113   {
2114     return aRetVal;
2115   }
2116
2117   // loop by edges
2118   TopTools_IndexedDataMapOfShapeListOfShape aMapEF;
2119   TopExp::MapShapesAndAncestors(theShape,TopAbs_EDGE,TopAbs_FACE,aMapEF);
2120   for(Standard_Integer anInd = 1; anInd <= aMapEF.Extent(); anInd++)
2121   {
2122     const TopoDS_Edge& anEdg = TopoDS::Edge(aMapEF.FindKey(anInd));
2123     const TopTools_ListOfShape& anEdgList = aMapEF.FindFromIndex(anInd);
2124     if (anEdgList.Extent() != 2)
2125       continue;
2126     TopTools_ListIteratorOfListOfShape anItF(anEdgList);
2127     const TopoDS_Face aFace1 = TopoDS::Face(anItF.Value());
2128     anItF.Next();
2129     const TopoDS_Face aFace2 = TopoDS::Face(anItF.Value());
2130     TopLoc_Location aLoc1, aLoc2;
2131     const Handle(Poly_Triangulation)& aPT1 = BRep_Tool::Triangulation(aFace1, aLoc1);
2132     const Handle(Poly_Triangulation)& aPT2 = BRep_Tool::Triangulation(aFace2, aLoc2);
2133
2134     if(aPT1.IsNull() || aPT2.IsNull())
2135       continue;
2136
2137     if(!aPT1->HasNormals() || !aPT2->HasNormals())
2138       continue;
2139
2140     const Handle(Poly_PolygonOnTriangulation)& aPTEF1 = 
2141                                 BRep_Tool::PolygonOnTriangulation(anEdg, aPT1, aLoc1);
2142     const Handle(Poly_PolygonOnTriangulation)& aPTEF2 = 
2143                                 BRep_Tool::PolygonOnTriangulation(anEdg, aPT2, aLoc2);
2144
2145     TShort_Array1OfShortReal& aNormArr1 = aPT1->ChangeNormals();
2146     TShort_Array1OfShortReal& aNormArr2 = aPT2->ChangeNormals();
2147
2148     if (aPTEF1->Nodes().Lower() != aPTEF2->Nodes().Lower() || 
2149         aPTEF1->Nodes().Upper() != aPTEF2->Nodes().Upper()) 
2150       continue; 
2151
2152     for(Standard_Integer anEdgNode = aPTEF1->Nodes().Lower();
2153                          anEdgNode <= aPTEF1->Nodes().Upper(); anEdgNode++)
2154     {
2155       //Number of node
2156       const Standard_Integer aFNodF1 = aPTEF1->Nodes().Value(anEdgNode);
2157       const Standard_Integer aFNodF2 = aPTEF2->Nodes().Value(anEdgNode);
2158
2159       const Standard_Integer aFNorm1FirstIndex = aNormArr1.Lower() + 3*
2160                                                     (aFNodF1 - aPT1->Nodes().Lower());
2161       const Standard_Integer aFNorm2FirstIndex = aNormArr2.Lower() + 3*
2162                                                     (aFNodF2 - aPT2->Nodes().Lower());
2163
2164       gp_XYZ aNorm1(aNormArr1.Value(aFNorm1FirstIndex),
2165                     aNormArr1.Value(aFNorm1FirstIndex+1),
2166                     aNormArr1.Value(aFNorm1FirstIndex+2));
2167       gp_XYZ aNorm2(aNormArr2.Value(aFNorm2FirstIndex),
2168                     aNormArr2.Value(aFNorm2FirstIndex+1),
2169                     aNormArr2.Value(aFNorm2FirstIndex+2));
2170       const Standard_Real aDot = aNorm1 * aNorm2;
2171
2172       if(aDot > aThresDot)
2173       {
2174         gp_XYZ aNewNorm = (aNorm1 + aNorm2).Normalized();
2175         aNormArr1.ChangeValue(aFNorm1FirstIndex) =
2176           aNormArr2.ChangeValue(aFNorm2FirstIndex) =
2177           static_cast<Standard_ShortReal>(aNewNorm.X());
2178         aNormArr1.ChangeValue(aFNorm1FirstIndex+1) =
2179           aNormArr2.ChangeValue(aFNorm2FirstIndex+1) =
2180           static_cast<Standard_ShortReal>(aNewNorm.Y());
2181         aNormArr1.ChangeValue(aFNorm1FirstIndex+2) =
2182           aNormArr2.ChangeValue(aFNorm2FirstIndex+2) =
2183           static_cast<Standard_ShortReal>(aNewNorm.Z());
2184         aRetVal = Standard_True;
2185       }
2186     }
2187   }
2188
2189   return aRetVal;
2190 }
2191
2192 //=======================================================================
2193 //function : SortFaces
2194 //purpose  : 
2195 //=======================================================================
2196
2197 void  BRepLib::SortFaces (const TopoDS_Shape& Sh,
2198   TopTools_ListOfShape& LF)
2199 {
2200   LF.Clear();
2201   TopTools_ListOfShape LTri,LPlan,LCyl,LCon,LSphere,LTor,LOther;
2202   TopExp_Explorer exp(Sh,TopAbs_FACE);
2203   TopLoc_Location l;
2204   Handle(Geom_Surface) S;
2205
2206   for (; exp.More(); exp.Next()) {
2207     const TopoDS_Face&   F = TopoDS::Face(exp.Current());
2208     S = BRep_Tool::Surface(F, l);
2209     if (!S.IsNull()) {
2210       if (S->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
2211         S = Handle(Geom_RectangularTrimmedSurface)::DownCast (S)->BasisSurface();
2212       }
2213       GeomAdaptor_Surface AS(S);
2214       switch (AS.GetType()) {
2215       case GeomAbs_Plane: 
2216         {
2217           LPlan.Append(F);
2218           break;
2219         }
2220       case GeomAbs_Cylinder: 
2221         {
2222           LCyl.Append(F);
2223           break;
2224         }
2225       case GeomAbs_Cone: 
2226         {
2227           LCon.Append(F);
2228           break;
2229         }
2230       case GeomAbs_Sphere: 
2231         {
2232           LSphere.Append(F);
2233           break;
2234         }
2235       case GeomAbs_Torus: 
2236         {
2237           LTor.Append(F);
2238           break;
2239         }
2240       default:
2241         LOther.Append(F);
2242       }
2243     }
2244     else LTri.Append(F);
2245   }
2246   LF.Append(LPlan); LF.Append(LCyl  ); LF.Append(LCon); LF.Append(LSphere);
2247   LF.Append(LTor ); LF.Append(LOther); LF.Append(LTri); 
2248 }
2249
2250 //=======================================================================
2251 //function : ReverseSortFaces
2252 //purpose  : 
2253 //=======================================================================
2254
2255 void  BRepLib::ReverseSortFaces (const TopoDS_Shape& Sh,
2256   TopTools_ListOfShape& LF)
2257 {
2258   LF.Clear();
2259   // Use the allocator of the result LF for intermediate results
2260   TopTools_ListOfShape LTri(LF.Allocator()), LPlan(LF.Allocator()),
2261     LCyl(LF.Allocator()), LCon(LF.Allocator()), LSphere(LF.Allocator()),
2262     LTor(LF.Allocator()), LOther(LF.Allocator());
2263   TopExp_Explorer exp(Sh,TopAbs_FACE);
2264   TopLoc_Location l;
2265
2266   for (; exp.More(); exp.Next()) {
2267     const TopoDS_Face&   F = TopoDS::Face(exp.Current());
2268     const Handle(Geom_Surface)& S = BRep_Tool::Surface(F, l);
2269     if (!S.IsNull()) {
2270       GeomAdaptor_Surface AS(S);
2271       switch (AS.GetType()) {
2272       case GeomAbs_Plane: 
2273         {
2274           LPlan.Append(F);
2275           break;
2276         }
2277       case GeomAbs_Cylinder: 
2278         {
2279           LCyl.Append(F);
2280           break;
2281         }
2282       case GeomAbs_Cone: 
2283         {
2284           LCon.Append(F);
2285           break;
2286         }
2287       case GeomAbs_Sphere: 
2288         {
2289           LSphere.Append(F);
2290           break;
2291         }
2292       case GeomAbs_Torus: 
2293         {
2294           LTor.Append(F);
2295           break;
2296         }
2297       default:
2298         LOther.Append(F);
2299       }
2300     }
2301     else LTri.Append(F);
2302   }
2303   LF.Append(LTri); LF.Append(LOther); LF.Append(LTor ); LF.Append(LSphere);
2304   LF.Append(LCon); LF.Append(LCyl  ); LF.Append(LPlan);
2305
2306 }
2307
2308 //=======================================================================
2309 // function: BoundingVertex
2310 // purpose : 
2311 //=======================================================================
2312 void BRepLib::BoundingVertex(const NCollection_List<TopoDS_Shape>& theLV,
2313                              gp_Pnt& theNewCenter, Standard_Real& theNewTol)
2314 {
2315   Standard_Integer aNb;
2316   //
2317   aNb=theLV.Extent();
2318   if (aNb < 2) {
2319     return;
2320   }
2321   //
2322   else if (aNb==2) {
2323     Standard_Integer m, n;
2324     Standard_Real aR[2], dR, aD, aEps;
2325     TopoDS_Vertex aV[2];
2326     gp_Pnt aP[2];
2327     //
2328     aEps=RealEpsilon();
2329     for (m=0; m<aNb; ++m) {
2330       aV[m]=(!m)? 
2331         *((TopoDS_Vertex*)(&theLV.First())):
2332         *((TopoDS_Vertex*)(&theLV.Last()));
2333       aP[m]=BRep_Tool::Pnt(aV[m]);
2334       aR[m]=BRep_Tool::Tolerance(aV[m]);
2335     }  
2336     //
2337     m=0; // max R
2338     n=1; // min R
2339     if (aR[0]<aR[1]) {
2340       m=1;
2341       n=0;
2342     }
2343     //
2344     dR=aR[m]-aR[n]; // dR >= 0.
2345     gp_Vec aVD(aP[m], aP[n]);
2346     aD=aVD.Magnitude();
2347     //
2348     if (aD<=dR || aD<aEps) { 
2349       theNewCenter = aP[m];
2350       theNewTol = aR[m];
2351     }
2352     else {
2353       Standard_Real aRr;
2354       gp_XYZ aXYZr;
2355       gp_Pnt aPr;
2356       //
2357       aRr=0.5*(aR[m]+aR[n]+aD);
2358       aXYZr=0.5*(aP[m].XYZ()+aP[n].XYZ()-aVD.XYZ()*(dR/aD));
2359       aPr.SetXYZ(aXYZr);
2360       //
2361       theNewCenter = aPr;
2362       theNewTol = aRr;
2363       //aBB.MakeVertex (aVnew, aPr, aRr);
2364     }
2365     return;
2366   }// else if (aNb==2) {
2367   //
2368   else { // if (aNb>2)
2369     // compute the point
2370     //
2371     // issue 0027540 - sum of doubles may depend on the order
2372     // of addition, thus sort the coordinates for stable result
2373     Standard_Integer i;
2374     NCollection_Array1<gp_Pnt> aPoints(0, aNb-1);
2375     NCollection_List<TopoDS_Edge>::Iterator aIt(theLV);
2376     for (i = 0; aIt.More(); aIt.Next(), ++i) {
2377       const TopoDS_Vertex& aVi = *((TopoDS_Vertex*)(&aIt.Value()));
2378       gp_Pnt aPi = BRep_Tool::Pnt(aVi);
2379       aPoints(i) = aPi;
2380     }
2381     //
2382     std::sort(aPoints.begin(), aPoints.end(), BRepLib_ComparePoints());
2383     //
2384     gp_XYZ aXYZ(0., 0., 0.);
2385     for (i = 0; i < aNb; ++i) {
2386       aXYZ += aPoints(i).XYZ();
2387     }
2388     aXYZ.Divide((Standard_Real)aNb);
2389     //
2390     gp_Pnt aP(aXYZ);
2391     //
2392     // compute the tolerance for the new vertex
2393     Standard_Real aTi, aDi, aDmax;
2394     //
2395     aDmax=-1.;
2396     aIt.Initialize(theLV);
2397     for (; aIt.More(); aIt.Next()) {
2398       TopoDS_Vertex& aVi=*((TopoDS_Vertex*)(&aIt.Value()));
2399       gp_Pnt aPi=BRep_Tool::Pnt(aVi);
2400       aTi=BRep_Tool::Tolerance(aVi);
2401       aDi=aP.SquareDistance(aPi);
2402       aDi=sqrt(aDi);
2403       aDi=aDi+aTi;
2404       if (aDi > aDmax) {
2405         aDmax=aDi;
2406       }
2407     }
2408     //
2409     theNewCenter = aP;
2410     theNewTol = aDmax;
2411   }
2412 }