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