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