0022885: Bugfix: else clause applies to the wrong if statement because of missing...
[occt.git] / src / ShapeAnalysis / ShapeAnalysis_CheckSmallFace.cxx
1 #include <ShapeAnalysis_CheckSmallFace.ixx>
2 #include <Standard_ErrorHandler.hxx>  
3 #include <TopTools_ListOfShape.hxx>
4 #include <TColStd_ListOfReal.hxx>
5 #include <ShapeExtend.hxx>
6 #include <gp_Pnt.hxx>
7 #include <TopoDS_Vertex.hxx>
8 #include <TopExp_Explorer.hxx>
9 #include <TopoDS.hxx>
10 #include <TopoDS_Shell.hxx>
11 #include <TopoDS_Edge.hxx>
12 #include <TopoDS_Vertex.hxx>
13 #include <TopoDS_Iterator.hxx>
14 #include <TopExp_Explorer.hxx>
15 #include <TopExp.hxx>
16 #include <BRep_Tool.hxx>
17
18 #include <Geom_Curve.hxx>
19 #include <Geom_Surface.hxx>
20 #include <Geom_BSplineSurface.hxx>
21 #include <Geom_BezierSurface.hxx>
22
23 #include <ShapeAnalysis_Curve.hxx>
24 //#include <GeomLProp_SLProps.hxx>
25 #include <GeomAdaptor_Surface.hxx>
26 #include <Geom_ElementarySurface.hxx>
27 #include <TColStd_Array2OfReal.hxx>
28 #include <gp_Dir.hxx>
29 #include <gp_Vec.hxx>
30 #include <TColStd_Array1OfReal.hxx>
31 #include <TColgp_Array1OfPnt.hxx>
32 #include <TopTools_Array1OfShape.hxx>
33
34 #include <TColgp_Array2OfPnt.hxx>
35
36 #include <TColgp_SequenceOfXYZ.hxx>
37 #include <BRep_Builder.hxx>
38 #include <BRepTools.hxx>
39 #include <Precision.hxx>
40 #include <TopoDS_Wire.hxx>
41 //#include <ShapeFix_Wire.hxx>
42 #include <Geom_Line.hxx>
43 #include <TopExp.hxx>
44 #include <ElCLib.hxx>
45 #include <TopoDS_Builder.hxx>
46 #include <TopoDS_Compound.hxx>
47 #include <Poly_Polygon3D.hxx>
48 #include <Geom2d_Curve.hxx>
49 #include <BRepLib.hxx>
50 #include <GeomLib.hxx>
51
52 #include <Geom_TrimmedCurve.hxx>
53 #include <GeomAdaptor_Curve.hxx>
54 #include <ShapeAnalysis_WireOrder.hxx>
55 #include <ShapeExtend_WireData.hxx>
56 #include <ShapeAnalysis_Wire.hxx>
57 #include <TopTools_HSequenceOfShape.hxx>
58 #include <TopoDS_Iterator.hxx> 
59 //=======================================================
60 //function : ShapeAnalysis_CheckSmallFace
61 //purpose  : 
62 //=======================================================================
63
64 ShapeAnalysis_CheckSmallFace::ShapeAnalysis_CheckSmallFace()
65 {
66   myStatusSpot = ShapeExtend::EncodeStatus ( ShapeExtend_OK );
67   myStatusStrip = ShapeExtend::EncodeStatus ( ShapeExtend_OK );
68   myStatusPin = ShapeExtend::EncodeStatus ( ShapeExtend_OK );
69   myStatusTwisted = ShapeExtend::EncodeStatus ( ShapeExtend_OK );
70   myStatusSplitVert = ShapeExtend::EncodeStatus ( ShapeExtend_OK );
71   
72 }
73 static void MinMaxPnt
74   (const gp_Pnt& p, Standard_Integer& nb,
75    Standard_Real& minx, Standard_Real& miny, Standard_Real& minz,
76    Standard_Real& maxx, Standard_Real& maxy, Standard_Real& maxz)
77 {
78   Standard_Real x,y,z;
79   p.Coord (x,y,z);
80   if (nb < 1)  {  minx = maxx = x; miny = maxy = y; minz = maxz = z;  }
81   else  {
82     if (minx > x) minx = x;  if (maxx < x) maxx = x;
83     if (miny > y) miny = y;  if (maxy < y) maxy = y;
84     if (minz > z) minz = z;  if (maxz < z) maxz = z;
85   }
86   nb ++;
87 }
88
89 static Standard_Boolean MinMaxSmall
90 (const Standard_Real minx, const Standard_Real miny, const Standard_Real minz, const Standard_Real maxx, const Standard_Real maxy, const Standard_Real maxz, const Standard_Real toler)
91 {
92   Standard_Real dx = maxx - minx;
93   Standard_Real dy = maxy - miny;
94   Standard_Real dz = maxz - minz;
95
96   if ((dx > toler && !Precision::IsInfinite (dx)) ||
97       (dy > toler && !Precision::IsInfinite (dy)) ||
98       (dz > toler && !Precision::IsInfinite (dz)))
99     return Standard_False;
100   return Standard_True;
101 }
102
103 //=======================================================================
104 //function : IsSpotFace
105 //purpose  : 
106 //=======================================================================
107
108  Standard_Integer ShapeAnalysis_CheckSmallFace::IsSpotFace(const TopoDS_Face& F,gp_Pnt& spot,Standard_Real& spotol,const Standard_Real tol) const
109 {
110
111   Standard_Real toler = tol;  Standard_Real tolv = tol;
112 //  Compute tolerance to get : from greatest tol of vertices
113 //  In addition, also computes min-max of vertices
114 //  To finally compare mini-max box with tolerance
115   // gka Mar2000 Protection against faces without wires
116   // but they occur due to bugs in the algorithm itself, it needs to be fixed
117   Standard_Boolean isWir = Standard_False;
118   for(TopoDS_Iterator itw(F,Standard_False) ; itw.More();itw.Next()) {
119     if(itw.Value().ShapeType() != TopAbs_WIRE)
120       continue;
121     TopoDS_Wire w1 = TopoDS::Wire(itw.Value());
122     if (!w1.IsNull()) {isWir = Standard_True; break;}
123   }
124   if(!isWir) return Standard_True;
125   Standard_Integer nbv = 0;
126   Standard_Real minx =0 ,miny = 0 ,minz = 0,maxx = Precision::Infinite(), maxy = Precision::Infinite(),maxz = Precision::Infinite();
127   TopoDS_Vertex V0;
128   Standard_Boolean same = Standard_True;
129   for (TopExp_Explorer iv(F,TopAbs_VERTEX); iv.More(); iv.Next()) {
130     TopoDS_Vertex V = TopoDS::Vertex (iv.Current());
131     if (V0.IsNull()) V0 = V;
132     else if (same) { if (!V0.IsSame(V)) same = Standard_False; }
133
134     gp_Pnt pnt = BRep_Tool::Pnt (V);
135     // Standard_Real x,y,z;
136     MinMaxPnt (pnt, nbv, minx,miny,minz, maxx,maxy,maxz);
137
138     if (tol < 0) {
139       tolv = BRep_Tool::Tolerance (V);
140       if (tolv > toler) toler = tolv;
141     }
142   }
143
144 //   Now, testing
145   if (!MinMaxSmall(minx,miny,minz,maxx,maxy,maxz,toler)) return 0;
146
147 //   All vertices are confused
148 //   Check edges (a closed edge may be a non-null length edge !)
149 //   By picking intermediate point on each one
150   for (TopExp_Explorer ie(F,TopAbs_EDGE); ie.More(); ie.Next()) {
151     TopoDS_Edge E = TopoDS::Edge (ie.Current());
152     Standard_Real cf,cl;
153     Handle(Geom_Curve) C3D = BRep_Tool::Curve (E,cf,cl);
154     if (C3D.IsNull()) continue;
155     gp_Pnt debut  = C3D->Value (cf);
156     gp_Pnt milieu = C3D->Value ( (cf+cl)/2);
157     if (debut.SquareDistance(milieu) > toler*toler) return 0;
158   }
159
160   spot.SetCoord ( (minx+maxx)/2. , (miny+maxy)/2. , (minz+maxz)/2. );
161   spotol = maxx-minx;
162   spotol = Max (spotol, maxy-miny);
163   spotol = Max (spotol, maxz-minz);
164   spotol = spotol/2.;
165
166   return (same ? 2 : 1);
167 }
168
169 //=======================================================================
170 //function : CheckSpotFace
171 //purpose  : 
172 //=======================================================================
173
174  Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckSpotFace(const TopoDS_Face& F,const Standard_Real tol) 
175 {
176   gp_Pnt spot;
177   Standard_Real spotol;
178   Standard_Integer stat = IsSpotFace (F,spot,spotol,tol);
179   if(!stat) return Standard_False; 
180   switch(stat) {
181     case 1: myStatusSpot = ShapeExtend::EncodeStatus (ShapeExtend_DONE1); break;
182     case 2: myStatusSpot = ShapeExtend::EncodeStatus (ShapeExtend_DONE2); break;
183     default : break;
184   
185   }
186   return Standard_True;
187 }
188
189 //=======================================================================
190 //function : IsStripSupport
191 //purpose  : 
192 //=======================================================================
193
194  Standard_Boolean ShapeAnalysis_CheckSmallFace::IsStripSupport(const TopoDS_Face& F,const Standard_Real tol) 
195 {
196
197   Standard_Real toler = tol;
198   if (toler < 0) toler = 1.e-07;    // ?? better to compute tolerance zones
199
200   TopLoc_Location loc;
201   Handle(Geom_Surface) surf = BRep_Tool::Surface (F,loc);
202   if (surf.IsNull()) return 0;
203
204 //  Checking on poles for bezier-bspline
205 //  A more general way is to check Values by scanning ISOS (slower)
206
207   Handle(Geom_BSplineSurface) bs = Handle(Geom_BSplineSurface)::DownCast(surf);
208   Handle(Geom_BezierSurface)  bz = Handle(Geom_BezierSurface)::DownCast(surf);
209
210   // Standard_Integer stat = 2;  // 2 : small in V direction
211   if (!bs.IsNull() || !bz.IsNull()) {
212     Standard_Boolean cbz = (!bz.IsNull());
213     Standard_Integer iu,iv, nbu, nbv;
214     if (cbz) { nbu = bz->NbUPoles(), nbv = bz->NbVPoles(); }
215     else     { nbu = bs->NbUPoles(), nbv = bs->NbVPoles(); }
216     // Standard_Real dx = 0, dy = 0, dz = 0;
217     // Standard_Real    x,y,z;
218     Standard_Real minx,miny,minz,maxx,maxy,maxz;
219     Standard_Boolean issmall = Standard_True;
220
221     for (iu = 1; iu <= nbu; iu ++) {
222 //    for each U line, scan poles in V (V direction)
223       Standard_Integer nb = 0;
224       for (iv = 1; iv <= nbv; iv ++) {
225         gp_Pnt unp = (cbz ? bz->Pole(iu,iv) : bs->Pole(iu,iv));
226         MinMaxPnt (unp, nb, minx,miny,minz, maxx,maxy,maxz);
227       }
228       if (!MinMaxSmall(minx,miny,minz,maxx,maxy,maxz,toler))
229         {  issmall = Standard_False;  break;  }    // small in V ?
230     }
231     if (issmall) {
232       myStatusStrip = ShapeExtend::EncodeStatus ( ShapeExtend_DONE2);
233       return issmall;    // OK, small in V
234     }
235     issmall = Standard_True;
236     for (iv = 1; iv <= nbv; iv ++) {
237 //    for each V line, scan poles in U (U direction)
238       Standard_Integer nb = 0;
239       for (iu = 1; iu <= nbu; iu ++) {
240         gp_Pnt unp = (cbz ? bz->Pole(iu,iv) : bs->Pole(iu,iv));
241         MinMaxPnt (unp, nb, minx,miny,minz, maxx,maxy,maxz);
242       }
243       if (!MinMaxSmall(minx,miny,minz,maxx,maxy,maxz,toler))
244         {  issmall = Standard_False;  break;  }    // small in U ?
245     }
246     if (issmall) {
247       myStatusStrip = ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
248       return issmall;
249     }// OK, small in U
250   }
251
252   return Standard_False;
253 }
254
255 //=======================================================================
256 //function : CheckStripEdges
257 //purpose  : 
258 //=======================================================================
259
260  Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckStripEdges(const TopoDS_Edge& E1,const TopoDS_Edge& E2,const Standard_Real tol,Standard_Real& dmax) const
261 {
262   //  We have the topological configuration OK : 2 edges, 2 vertices
263   //  But, are these two edges well confused ?
264   Standard_Real toler = tol;
265   if (tol < 0) {
266     Standard_Real tole = BRep_Tool::Tolerance(E1) + BRep_Tool::Tolerance(E2);
267     if (toler < tole / 2.) toler = tole/2.;
268   }
269   
270   //   We project a list of points from each curve, on the opposite one,
271   //   we check the distance
272   Standard_Integer nbint = 10;
273   
274   ShapeAnalysis_Curve SAC;
275   Standard_Real cf1,cl1,cf2,cl2,u; dmax = 0;
276   Handle(Geom_Curve) C1,C2;
277   C1 = BRep_Tool::Curve (E1,cf1,cl1);
278   C2 = BRep_Tool::Curve (E2,cf2,cl2);
279   if(C1.IsNull() || C2.IsNull()) return Standard_False;
280   cf1 = Max(cf1, C1->FirstParameter());
281   cl1 = Min(cl1, C1->LastParameter());
282   Handle(Geom_TrimmedCurve) C1T = new Geom_TrimmedCurve(C1,cf1,cl1,Standard_True);
283   //pdn protection against feature in Trimmed_Curve
284   cf1 = C1T->FirstParameter();
285   cl1 = C1T->LastParameter();
286   Handle(Geom_TrimmedCurve) CC;
287   cf2 = Max(cf2, C2->FirstParameter());
288   cl2 = Min(cl2, C2->LastParameter());
289   Handle(Geom_TrimmedCurve) C2T = new Geom_TrimmedCurve(C2,cf2,cl2, Standard_True);
290   cf2 = C2T->FirstParameter();
291   cl2 = C2T->LastParameter();
292   
293   Standard_Real cd1 = (cl1 - cf1)/nbint;
294   Standard_Real cd2 = (cl2 - cf2)/nbint;
295   Standard_Real f,l;
296   f = cf2; l = cl2;
297   for (int numcur = 0; numcur < 2; numcur ++) {
298     u = cf1;
299     if (numcur)  {  CC = C1T; C1T = C2T; C2T = CC;
300                     cd1 = cd2; //smh added replacing step and replacing first
301                     u = cf2;   //parameter
302                     f = cf1; l = cl1;
303                   }
304     for (int nump = 0; nump <= nbint; nump ++) {
305       gp_Pnt p2, p1 = C1T->Value (u);
306       Standard_Real para;
307       //pdn Adaptor curve is used to avoid of enhancing of domain.
308       GeomAdaptor_Curve GAC(C2T);
309       Standard_Real dist = SAC.Project (GAC,p1,toler,p2,para);
310       //pdn check if parameter of projection is in the domain of the edge.
311       if (para < f || para > l) return Standard_False;
312       if (dist > dmax) dmax = dist;
313       if (dist > toler) return Standard_False;
314       u += cd1;
315     }
316   }
317   return (dmax < toler);
318 }
319
320 //=======================================================================
321 //function : FindStripEdges
322 //purpose  : 
323 //=======================================================================
324
325  Standard_Boolean ShapeAnalysis_CheckSmallFace::FindStripEdges(const TopoDS_Face& F,TopoDS_Edge& E1,TopoDS_Edge& E2,const Standard_Real tol,Standard_Real& dmax) 
326 {
327   E1.Nullify();  E2.Nullify();
328   Standard_Integer nb = 0;
329   for (TopExp_Explorer ex(F,TopAbs_EDGE); ex.More(); ex.Next()) {
330     TopoDS_Edge E = TopoDS::Edge (ex.Current());
331     if (nb == 1 && E.IsSame(E1))
332       continue; // ignore seam edge
333     TopoDS_Vertex V1,V2;
334     TopExp::Vertices (E,V1,V2);
335     gp_Pnt p1,p2;
336     p1 = BRep_Tool::Pnt (V1);
337     p2 = BRep_Tool::Pnt (V2);
338     Standard_Real toler = tol;
339     if (toler <= 0) toler = (BRep_Tool::Tolerance(V1) + BRep_Tool::Tolerance(V2) ) / 2.;
340
341 //    Extremities
342     Standard_Real dist = p1.Distance(p2);
343 //    Middle point
344     Standard_Real cf,cl;
345     Handle(Geom_Curve) CC;
346     CC = BRep_Tool::Curve (E,cf,cl);
347     Standard_Boolean isNullLength = Standard_True;
348     if (!CC.IsNull()) {
349       gp_Pnt pp = CC->Value ( (cf+cl)/2.);
350       if (pp.Distance(p1) < toler && pp.Distance(p2) < toler) continue;
351       isNullLength = Standard_False;
352     }
353     if (dist <= toler && isNullLength) continue; //smh
354     nb ++;
355     if (nb == 1) E1 = E;
356     else if (nb == 2) E2 = E;
357     else return Standard_False;
358   }
359   //   Now, check these two edge to define a strip !
360   if (!E1.IsNull()&&!E2.IsNull()) 
361     if(!CheckStripEdges (E1,E2,tol,dmax)) return Standard_False; 
362       else {   
363         myStatusStrip = ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
364         return Standard_True ;
365       }
366   return Standard_False;
367 }
368
369 //=======================================================================
370 //function : CheckSingleStrip
371 //purpose  : 
372 //=======================================================================
373
374  Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckSingleStrip(const TopoDS_Face& F,
375                                                                 TopoDS_Edge& E1, TopoDS_Edge& E2,const Standard_Real tol) 
376 {
377  Standard_Real toler = tol;
378   Standard_Real minx,miny,minz,maxx,maxy,maxz;
379
380 // In this case, we have 2 vertices and 2 great edges. Plus possibly 2 small
381 //    edges, one on each vertex
382   TopoDS_Vertex V1,V2;
383   Standard_Integer nb = 0;
384   for (TopExp_Explorer itv (F,TopAbs_VERTEX); itv.More(); itv.Next()) {
385     TopoDS_Vertex V = TopoDS::Vertex (itv.Current());
386     if (V1.IsNull()) V1 = V;
387     else if (V1.IsSame(V)) continue;
388     else if (V2.IsNull()) V2 = V;
389     else if (V2.IsSame(V)) continue;
390     else return 0;
391   }
392
393 // Checking edges
394   //TopoDS_Edge E1,E2;
395   nb = 0;
396   for (TopExp_Explorer ite (F,TopAbs_EDGE); ite.More(); ite.Next()) {
397     TopoDS_Edge E = TopoDS::Edge (ite.Current());
398     if (nb == 1 && E.IsSame(E1))
399       continue; // ignore seam edge
400     TopoDS_Vertex VA,VB;
401     TopExp::Vertices (E,VA,VB);
402     if (tol < 0) {
403       Standard_Real tolv;
404       tolv = BRep_Tool::Tolerance (VA);
405       if (toler < tolv) toler = tolv;
406       tolv = BRep_Tool::Tolerance (VB);
407       if (toler < tolv) toler = tolv;
408     }
409
410 //    Edge on same vertex : small one ?
411     if (VA.IsSame(VB)) {
412       Standard_Real cf,cl;
413       Handle(Geom_Curve) C3D;
414       if (!BRep_Tool::Degenerated(E)) C3D = BRep_Tool::Curve (E,cf,cl);
415       if (C3D.IsNull()) continue;  // DGNR
416       Standard_Integer np = 0;
417       gp_Pnt deb = C3D->Value(cf);
418       MinMaxPnt (deb,np,minx,miny,minz,maxx,maxy,maxz);
419       gp_Pnt fin = C3D->Value(cl);
420       MinMaxPnt (fin,np,minx,miny,minz,maxx,maxy,maxz);
421       gp_Pnt mid = C3D->Value( (cf+cl)/2. );
422       MinMaxPnt (mid,np,minx,miny,minz,maxx,maxy,maxz);
423       if (!MinMaxSmall (minx,miny,minz,maxx,maxy,maxz,toler)) return Standard_False;
424     } else {
425 //    Other case : two maximum allowed
426       nb ++;
427       if (nb > 2) return Standard_False;
428       if (nb == 1)  {  V1 = VA;  V2 = VB;  E1 = E;  }
429       else if (nb == 2) {
430         if (V1.IsSame(VA) && !V2.IsSame(VB)) return Standard_False;
431         if (V1.IsSame(VB) && !V2.IsSame(VA)) return Standard_False;
432         E2 = E;
433       }
434       else return Standard_False;
435     }
436   }
437
438   if (nb < 2) return Standard_False;   // only one vertex : cannot be a strip ...
439
440 //   Checking if E1 and E2 define a Strip
441   Standard_Real dmax;
442  if (!CheckStripEdges (E1,E2,tol,dmax)) return Standard_False;
443  myStatusStrip = ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
444  return Standard_True;
445 }
446
447 //=======================================================================
448 //function : CheckStripFace
449 //purpose  : 
450 //=======================================================================
451
452 Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckStripFace(const TopoDS_Face& F,
453                                                                TopoDS_Edge& E1, TopoDS_Edge& E2,const Standard_Real tol) 
454 {
455
456   // Standard_Integer stat;
457   if(CheckSingleStrip (F,E1,E2,tol)) return Standard_True ;   // it is a strip
458
459 //    IsStripSupport used as rejection. But this kind of test may be done
460 //    on ANY face, once we are SURE that FindStripEdges is reliable (and fast
461 //    enough)
462
463 //  ?? record a diagnostic StripFace, but without yet lists of edges
464 //  ??  Record Diagnostic "StripFace", no data (should be "Edges1" "Edges2")
465 //      but direction is known (1:U  2:V)
466    // TopoDS_Edge E1,E2;
467     Standard_Real dmax;
468     if(FindStripEdges (F,E1,E2,tol,dmax)) return Standard_True;
469
470 //  Now, trying edges : if there are 2 and only 2 edges greater than tolerance
471 //   (given or sum of vertex tolerances), do they define a strip
472 //  Warning : if yes, they bring different vertices ...
473
474   return Standard_False;
475
476 }
477
478 //=======================================================================
479 //function : CheckSplittingVertices
480 //purpose  : 
481 //=======================================================================
482
483  Standard_Integer ShapeAnalysis_CheckSmallFace::CheckSplittingVertices(const TopoDS_Face& F,
484                                                                        TopTools_DataMapOfShapeListOfShape& MapEdges,
485                                                                        ShapeAnalysis_DataMapOfShapeListOfReal& MapParam,
486                                                                        TopoDS_Compound& theAllVert)
487 {
488
489 //  Prepare array of vertices with their locations //TopTools
490   Standard_Integer nbv = 0, nbp = 0;
491   //TopoDS_Compound theAllVert;
492   BRep_Builder theBuilder;
493   //theBuilder.MakeCompound(theAllVert);
494   TopExp_Explorer itv; // svv Jan11 2000 : porting on DEC
495   for (itv.Init(F,TopAbs_VERTEX); itv.More(); itv.Next()) nbv ++;
496
497   if (nbv == 0) return 0;
498   TopTools_Array1OfShape vtx (1,nbv);
499   TColgp_Array1OfPnt vtp (1,nbv);
500   TColStd_Array1OfReal vto (1,nbv);
501
502   nbp = 0;
503   for (itv.Init(F,TopAbs_VERTEX); itv.More(); itv.Next()) {
504     nbp ++;
505     TopoDS_Vertex unv = TopoDS::Vertex (itv.Current());
506     vtx.SetValue (nbp,unv);
507     gp_Pnt unp = BRep_Tool::Pnt (unv);
508     vtp.SetValue (nbp,unp);
509     Standard_Real unt = myPrecision;
510     if (unt < 0) unt =BRep_Tool::Tolerance (unv);
511     vto.SetValue (nbp,unt);
512   }
513   nbv = nbp;  nbp = 0;  // now, counting splitting vertices
514
515 //  Check edges : are vertices (other than extremities) confused with it ?
516   ShapeAnalysis_Curve SAC;
517   for (Standard_Integer iv = 1; iv <= nbv; iv ++) {
518     TopoDS_Vertex V = TopoDS::Vertex (vtx.Value(iv));
519     TopTools_ListOfShape listEdge;
520     TColStd_ListOfReal listParam;
521     Standard_Boolean issplit = Standard_False;
522     for (TopExp_Explorer ite(F,TopAbs_EDGE); ite.More(); ite.Next()) {
523       TopoDS_Edge E = TopoDS::Edge (ite.Current());
524       TopoDS_Vertex V1,V2;
525       TopExp::Vertices (E,V1,V2);
526       Standard_Real cf,cl;
527       Handle(Geom_Curve) C3D = BRep_Tool::Curve (E,cf,cl);
528       if (C3D.IsNull()) continue;
529       if (V.IsSame(V1) || V.IsSame(V2)) continue;
530       gp_Pnt unp = vtp.Value(iv);
531       Standard_Real unt = vto.Value(iv);
532       gp_Pnt proj;
533       Standard_Real param;
534       Standard_Real dist = SAC.Project (C3D,unp,unt*10.,proj,param,cf,cl);
535       if (dist == 0.0) continue; //smh
536 //  Splitting Vertex to record ?
537       if (dist < unt) {
538 //  If Split occurs at beginning or end, it is not a split ...
539         Standard_Real fpar, lpar, eps = 1.e-06;
540         if (param >=cl || param <= cf) continue; // Out of range
541         fpar = param - cf; lpar = param - cl;
542         if ((Abs(fpar) < eps) || (Abs(lpar) < eps)) continue; // Near end or start 
543         listEdge.Append(E);
544         listParam.Append(param);
545         issplit = Standard_True;
546         
547       }
548     }
549     if(issplit) {
550       nbp ++;
551       theBuilder.Add(theAllVert, V);
552       MapEdges.Bind(V,listEdge);
553       MapParam.Bind(V,listParam);
554     }
555   }
556   if(nbp != 0) 
557   myStatusSplitVert = ShapeExtend::EncodeStatus (ShapeExtend_DONE);
558   return nbp;
559 }
560
561  
562 static Standard_Integer IsoStat
563   (const TColgp_Array2OfPnt& poles,
564    const Standard_Integer uorv, const Standard_Integer rank,
565    const Standard_Real tolpin, const Standard_Real toler)
566 {
567   Standard_Integer i, np = 0;
568   Standard_Integer i0 = (uorv == 1 ? poles.LowerCol() : poles.LowerRow());
569   Standard_Integer i1 = (uorv == 1 ? poles.UpperCol() : poles.UpperRow());
570   Standard_Real xmin,ymin,zmin, xmax,ymax,zmax;
571   for (i = i0; i <= i1; i ++) {
572     if (uorv == 1) MinMaxPnt (poles(rank,i),np,xmin,ymin,zmin, xmax,ymax,zmax);
573     else      MinMaxPnt (poles(i,rank), np, xmin,ymin,zmin, xmax,ymax,zmax);
574   }
575   if (MinMaxSmall (xmin,ymin,zmin, xmax,ymax,zmax, tolpin)) return 0;
576   if (MinMaxSmall (xmin,ymin,zmin, xmax,ymax,zmax, toler))  return 1;
577   return 2;
578 }
579
580 static Standard_Boolean CheckPoles(const TColgp_Array2OfPnt& poles, Standard_Integer uorv, Standard_Integer rank)
581 {
582   Standard_Integer i0 = (uorv == 1 ? poles.LowerCol() : poles.LowerRow());
583   Standard_Integer i1 = (uorv == 1 ? poles.UpperCol() : poles.UpperRow());
584   for (Standard_Integer i = i0; i <= i1-1; i ++) {
585     if (uorv == 1) {
586       if(poles(rank,i).IsEqual(poles(rank, i+1), 1e-15)) return Standard_True;
587     } else
588       if(poles(i,rank).IsEqual(poles(i+1,rank), 1e-15)) return Standard_True;
589   }  
590   return Standard_False;
591 }
592 //=======================================================================
593 //function : CheckPin
594 //purpose  : 
595 //=======================================================================
596
597 Standard_Boolean  ShapeAnalysis_CheckSmallFace::CheckPin (const TopoDS_Face& F, Standard_Integer& whatrow,Standard_Integer& sens)
598 {
599   TopLoc_Location loc;
600   Handle(Geom_Surface) surf = BRep_Tool::Surface (F,loc);
601   if (surf->IsKind(STANDARD_TYPE(Geom_ElementarySurface))) return Standard_False;
602
603   Standard_Real toler = myPrecision;
604   if (toler < 0) toler = 1.e-4;
605   Standard_Real tolpin = 1.e-9;  // for sharp sharp pin
606
607 //  Checking the poles
608
609 //  Take the poles : they give good idea of sharpness of a pin
610   Standard_Integer nbu = 0 , nbv = 0;
611   Handle(Geom_BSplineSurface) bs = Handle(Geom_BSplineSurface)::DownCast(surf);
612   Handle(Geom_BezierSurface)  bz = Handle(Geom_BezierSurface)::DownCast(surf);
613   if (!bs.IsNull())  {  nbu = bs->NbUPoles();  nbv = bs->NbVPoles();  }
614   if (!bz.IsNull())  {  nbu = bz->NbUPoles();  nbv = bz->NbVPoles();  }
615   if (nbu == 0 || nbv == 0) return Standard_False;
616
617   TColgp_Array2OfPnt allpoles (1,nbu,1,nbv);
618   if (!bs.IsNull()) bs->Poles (allpoles);
619   if (!bz.IsNull()) bz->Poles (allpoles);
620
621 //  Check each natural bound if it is a singularity (i.e. a pin)
622
623   sens = 0; 
624   Standard_Integer stat = 0;    // 0 none, 1 in U, 2 in V
625   whatrow = 0;  // 0 no row, else rank of row
626   stat = IsoStat(allpoles,1,  1,tolpin,toler);
627   if (stat) { sens = 1;  whatrow = nbu; }
628   
629   stat = IsoStat(allpoles,1,nbu,tolpin,toler);
630   if (stat) { sens = 1; whatrow = nbu;  }
631   
632   stat = IsoStat(allpoles,2,  1,tolpin,toler);
633   if (stat) { sens = 2; whatrow = 1; }
634   
635   stat = IsoStat(allpoles,2,nbv,tolpin,toler);
636   if (stat) { sens = 2; whatrow = nbv;  }
637
638   if (!sens) return Standard_False;    // no pin
639   
640   switch(stat) {
641   case 1: myStatusPin = ShapeExtend::EncodeStatus (ShapeExtend_DONE1); break;
642   case 2: myStatusPin = ShapeExtend::EncodeStatus (ShapeExtend_DONE2); break;
643     default : break;
644   }
645   //  cout<<(whatstat == 1 ? "Smooth" : "Sharp")<<" Pin on "<<(sens == 1 ? "U" : "V")<<" Row n0 "<<whatrow<<endl;
646   if (stat == 1 ) 
647     {
648       // Standard_Boolean EqualPoles = Standard_False;
649       if(CheckPoles(allpoles, 2, nbv)|| CheckPoles(allpoles, 2, 1)
650          ||CheckPoles(allpoles, 1, nbu)|| CheckPoles(allpoles, 1, 1))       
651         myStatusPin = ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
652     }
653   
654   
655   return Standard_True;
656 }
657
658 static Standard_Real TwistedNorm
659 (const Standard_Real x1, const Standard_Real y1, const Standard_Real z1, const Standard_Real x2, const Standard_Real y2, const Standard_Real z2)
660 {  return (x1*x2) + (y1*y2) + (z1*z2);  }
661
662 //=======================================================================
663 //function : CheckTwisted
664 //purpose  : 
665 //=======================================================================
666
667 Standard_Boolean  ShapeAnalysis_CheckSmallFace::CheckTwisted (const TopoDS_Face& F, Standard_Real& paramu,
668                                                              Standard_Real& paramv)
669 {
670   TopLoc_Location loc;
671   Handle(Geom_Surface) surf = BRep_Tool::Surface (F,loc);
672   if (surf->IsKind(STANDARD_TYPE(Geom_ElementarySurface))) return Standard_False;
673
674   Standard_Real toler = myPrecision;
675   if (toler < 0) toler = 1.e-4;
676 ////  GeomLProp_SLProps GLS (surf,2,toler);
677   GeomAdaptor_Surface GAS (surf);
678
679 // to be done : on isos of the surface
680 //  and on edges, at least of outer wire
681   Standard_Integer nbint = 5;
682   TColStd_Array2OfReal nx (1,nbint+1,1,nbint+1);
683   TColStd_Array2OfReal ny (1,nbint+1,1,nbint+1);
684   TColStd_Array2OfReal nz (1,nbint+1,1,nbint+1);
685   Standard_Integer iu,iv;
686   Standard_Real umin,umax,vmin,vmax;
687   surf->Bounds (umin,umax,vmin,vmax);
688   Standard_Real u = umin, du = (umax-umin)/nbint;
689   Standard_Real v = vmin, dv = (umax-umin)/nbint;
690
691   // gp_Dir norm;
692   for (iu = 1; iu <= nbint; iu ++) {
693     for (iv = 1; iv <= nbint; iv ++) {
694 //      GLS.SetParameters (u,v);
695 //      if (GLS.IsNormalDefined()) norm = GLS.Normal();
696       gp_Pnt curp;  gp_Vec V1,V2,VXnorm;
697       GAS.D1 (u,v,curp,V1,V2);
698       VXnorm = V1.Crossed(V2);
699       nx.SetValue (iu,iv,VXnorm.X());
700       ny.SetValue (iu,iv,VXnorm.Y());
701       nz.SetValue (iu,iv,VXnorm.Z());
702       v += dv;
703     }
704     u += du;
705     v = vmin;
706   }
707
708 //  Now, comparing normals on support surface, in both senses
709 //  In principle, it suffuces to check within outer bound
710
711   for (iu = 1; iu < nbint; iu ++) {
712     for (iv = 1; iv < nbint; iv ++) {
713 // We here check each normal (iu,iv) with (iu,iv+1) and with (iu+1,iv)
714 // if for each test, we have negative scalar product, this means angle > 90deg
715 // it is the criterion to say it is twisted
716       if (TwistedNorm ( nx(iu,iv),ny(iu,iv),nz(iu,iv) , nx(iu,iv+1),ny(iu,iv+1),nz(iu,iv+1) ) < 0. ||
717           TwistedNorm ( nx(iu,iv),ny(iu,iv),nz(iu,iv) , nx(iu+1,iv),ny(iu+1,iv),nz(iu+1,iv) ) < 0. ) {
718         myStatusTwisted = ShapeExtend::EncodeStatus (ShapeExtend_DONE);
719         paramu = umin+du*iu-du/2;
720         paramv = vmin+dv*iv-dv/2;
721         return Standard_True;
722       }
723     }
724   }
725
726 //   Now, comparing normals on edges ... to be done
727
728   return Standard_False;
729 }
730
731
732 //=======================================================================
733 //function : CheckPinFace
734 //purpose  : 
735 //=======================================================================
736 // Warning: This function not tested on many examples
737
738  Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckPinFace(const TopoDS_Face& F,
739   TopTools_DataMapOfShapeShape& mapEdges,const Standard_Real toler) 
740 {
741   //ShapeFix_Wire sfw;
742   TopExp_Explorer exp_w (F,TopAbs_WIRE); 
743   exp_w.More();
744   Standard_Real coef1=0, coef2; // =0 for deleting warning (skl)
745   TopoDS_Wire theCurWire = TopoDS::Wire (exp_w.Current());
746   ShapeAnalysis_WireOrder wi;
747   ShapeAnalysis_Wire sfw;
748   Handle(ShapeExtend_WireData) sbwd = new ShapeExtend_WireData(theCurWire);
749   sfw.Load(sbwd);
750   sfw.CheckOrder(wi);
751   Handle(TopTools_HSequenceOfShape) newedges = new TopTools_HSequenceOfShape();
752   Standard_Integer nb = wi.NbEdges();
753   Standard_Integer i = 0;
754   for ( i=1; i <= nb; i++ )
755     newedges->Append ( sbwd->Edge ( wi.Ordered(i) ) );
756   for ( i=1; i <= nb; i++ ) 
757     sbwd->Set ( TopoDS::Edge ( newedges->Value(i) ), i );
758   //sfw.Init(theCurWire,  F, Precision::Confusion());
759   //sfw.FixReorder();
760   //theCurWire = sfw.Wire();
761   theCurWire = sbwd->Wire();
762   i=1;
763   Standard_Boolean done = Standard_False;
764   Standard_Real tol = Precision::Confusion();
765   TopoDS_Edge theFirstEdge, theSecondEdge;
766   Standard_Real d1=0,d2=0;
767   for (TopExp_Explorer exp_e (F,TopAbs_EDGE); exp_e.More(); exp_e.Next()) 
768     {
769       TopoDS_Vertex V1,V2;
770       gp_Pnt p1, p2;
771       if (i==1) 
772         { 
773           theFirstEdge = TopoDS::Edge (exp_e.Current()); 
774           V1 = TopExp::FirstVertex(theFirstEdge);
775           V2 = TopExp::LastVertex(theFirstEdge);
776           p1 = BRep_Tool::Pnt(V1);
777           p2 = BRep_Tool::Pnt(V2);
778           tol = Max(BRep_Tool::Tolerance(V1), BRep_Tool::Tolerance(V2));
779           if (toler > 0) //tol = Max(tol, toler); gka
780           tol = toler;
781           d1 = p1.Distance(p2);
782           if (d1 == 0) return Standard_False;
783           if (d1/tol>=1) coef1 = d1/tol; else continue;
784           if (coef1<=3) continue;
785           i++; 
786           continue;
787         }
788       //Check the length of edge
789       theSecondEdge = TopoDS::Edge (exp_e.Current());
790       V1 = TopExp::FirstVertex(theSecondEdge);
791       V2 = TopExp::LastVertex(theSecondEdge);
792       
793       p1 = BRep_Tool::Pnt(V1);
794       p2 = BRep_Tool::Pnt(V2);
795       if (toler == -1) tol = Max(BRep_Tool::Tolerance(V1), BRep_Tool::Tolerance(V2)); 
796       else tol= toler;
797       if (p1.Distance(p2)> tol) continue;
798       //If there are two pin edges, record them in diagnostic
799       d2 = p1.Distance(p2); //gka
800       if (d2 == 0) return Standard_False;
801       if (d2/tol >= 1) coef2 = d2/tol; else continue;
802       if (coef2<=3) continue;
803       if (coef1>coef2*10) continue;
804       if (coef2>coef1*10)
805         {
806           theFirstEdge = theSecondEdge;
807           coef1 = coef2;
808           continue;
809         }
810       
811       if (CheckPinEdges(theFirstEdge, theSecondEdge, coef1, coef2,toler))
812         {
813           mapEdges.Bind(theFirstEdge,theSecondEdge);
814           myStatusPinFace = ShapeExtend::EncodeStatus (ShapeExtend_DONE);
815           done = Standard_True;
816         }     
817       
818       theFirstEdge = theSecondEdge;
819       coef1 = coef2;
820       //d1 = d2;
821     }
822   return done;
823 }
824
825
826 //=======================================================================
827 //function : CheckPinEdges
828 //purpose  : 
829 //=======================================================================
830 // Warning: This function not tested on many examples
831
832  Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckPinEdges(const TopoDS_Edge& theFirstEdge,const TopoDS_Edge& theSecondEdge,const Standard_Real coef1,
833    const Standard_Real coef2,const Standard_Real toler) const
834 {
835
836  Standard_Real cf1,cl1,cf2,cl2;
837  Handle(Geom_Curve) C1,C2,C3;
838  C1 = BRep_Tool::Curve (theFirstEdge,cf1,cl1);
839  C2 = BRep_Tool::Curve (theSecondEdge,cf2,cl2);
840  gp_Pnt p1, p2, pp1, pp2, pv;
841  Standard_Real d1 = (cf1-cl1)/coef1;
842  Standard_Real d2 = (cf2-cl2)/coef2;
843  //Standard_Real d1 = cf1-cl1/30; //10; gka
844   //Standard_Real d2 = cf2-cl2/30; //10;
845   p1 = C1->Value(cf1);
846   p2 = C1->Value(cl1);
847   pp1 = C2->Value(cf2);
848   pp2 = C2->Value(cl2);
849   Standard_Real tol;
850   Standard_Real paramc1=0, paramc2=0; // =0 for deleting warning (skl)
851   TopoDS_Vertex theSharedV = TopExp::LastVertex(theFirstEdge);
852   if  (toler == -1)  tol = BRep_Tool::Tolerance(theSharedV); else tol = toler;
853   pv = BRep_Tool::Pnt(theSharedV);
854   if (pv.Distance(p1)<=tol) paramc1 = cf1;
855   else if(pv.Distance(p2)<=tol) paramc1 = cl1;
856   if (pv.Distance(pp1)<=tol) paramc2 = cf2;
857   else if(pv.Distance(pp2)<=tol) paramc2 = cl2;
858   //Computing first derivative vectors and compare angle
859 //   gp_Vec V11, V12, V21, V22;
860 //   gp_Pnt tmp;
861 //   C1->D2(paramc1, tmp, V11, V21);
862 //   C2->D2(paramc2, tmp, V12, V22);
863 //   Standard_Real angle1, angle2;
864 //   try{
865 //     angle1 = V11.Angle(V12);
866 //     angle2 = V21.Angle(V22);
867 //   }
868 //   catch (Standard_Failure)
869 //     {
870 //       cout << "Couldn't compute angle between derivative vectors"  <<endl;
871 //       return Standard_False;
872 //     }
873 //   cout << "angle1 "   << angle1<< endl;
874 //   cout << "angle2 "   << angle2<< endl;
875 //   if (angle1<=0.0001) return Standard_True;
876   gp_Pnt proj;
877   if (p1.Distance(p2)<pp1.Distance(pp2)) 
878     {
879       C3=C1;
880       if (paramc1==cf1)
881        proj = C1->Value(paramc1 + (coef1-3)*d1);
882       else proj = C1->Value(paramc1-3*d1);
883         //proj = C1->Value(paramc1 + 9*d1);
884       //else proj = C1->Value(paramc1-d1);
885     }
886   else 
887     {
888       C3=C2;
889       if (paramc2==cf2)
890         proj = C2->Value(paramc2 + (coef2-3)*d2); 
891       else proj = C2->Value(paramc2 -3*d2);
892         //proj = C2->Value(paramc2 + 9*d2); 
893       //else proj = C2->Value(paramc2 -d2); 
894     }
895   Standard_Real param;
896   GeomAdaptor_Curve GAC(C3);
897   Standard_Real f = C3->FirstParameter();
898   Standard_Real l = C3->LastParameter();
899   gp_Pnt result;
900   ShapeAnalysis_Curve SAC;
901   Standard_Real dist = SAC.Project (GAC,proj,tol,result,param);
902   //pdn check if parameter of projection is in the domain of the edge.
903   if (param < f || param > l) return Standard_False;
904   if (dist > tol) return Standard_False;
905   if (dist <= tol) {
906      //Computing first derivative vectors and compare angle
907       gp_Vec V11, V12, V21, V22;
908       gp_Pnt tmp;
909       C1->D2(paramc1, tmp, V11, V21);
910       C2->D2(paramc2, tmp, V12, V22);
911       Standard_Real angle1=0, angle2=0;
912       try{
913         angle1 = V11.Angle(V12);
914         angle2 = V21.Angle(V22);
915       }
916       catch (Standard_Failure)
917         {
918           cout << "Couldn't compute angle between derivative vectors"  <<endl;
919           return Standard_False;
920         }
921 //       cout << "angle1 "   << angle1<< endl;
922 //       cout << "angle2 "   << angle2<< endl;
923       if ((angle1<=0.001 && angle2<=0.01) || ((M_PI-angle2)<= 0.001 && (M_PI-angle2)<= 0.01)) return Standard_True;
924       else return Standard_False;
925     } 
926     
927   return Standard_False;
928 }
929