e7fbcf899bb2ce3a0fd18acf7340dea6c86cd46d
[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   }
384   return Standard_False;
385 }
386
387 //=======================================================================
388 //function : CheckSingleStrip
389 //purpose  : 
390 //=======================================================================
391
392  Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckSingleStrip(const TopoDS_Face& F,
393                                                                 TopoDS_Edge& E1, TopoDS_Edge& E2,const Standard_Real tol) 
394 {
395  Standard_Real toler = tol;
396   Standard_Real minx,miny,minz,maxx,maxy,maxz;
397
398 // In this case, we have 2 vertices and 2 great edges. Plus possibly 2 small
399 //    edges, one on each vertex
400   TopoDS_Vertex V1,V2;
401   Standard_Integer nb = 0;
402   for (TopExp_Explorer itv (F,TopAbs_VERTEX); itv.More(); itv.Next()) {
403     TopoDS_Vertex V = TopoDS::Vertex (itv.Current());
404     if (V1.IsNull()) V1 = V;
405     else if (V1.IsSame(V)) continue;
406     else if (V2.IsNull()) V2 = V;
407     else if (V2.IsSame(V)) continue;
408     else return 0;
409   }
410
411 // Checking edges
412   //TopoDS_Edge E1,E2;
413   nb = 0;
414   for (TopExp_Explorer ite (F,TopAbs_EDGE); ite.More(); ite.Next()) {
415     TopoDS_Edge E = TopoDS::Edge (ite.Current());
416     if (nb == 1 && E.IsSame(E1))
417       continue; // ignore seam edge
418     TopoDS_Vertex VA,VB;
419     TopExp::Vertices (E,VA,VB);
420     if (tol < 0) {
421       Standard_Real tolv;
422       tolv = BRep_Tool::Tolerance (VA);
423       if (toler < tolv) toler = tolv;
424       tolv = BRep_Tool::Tolerance (VB);
425       if (toler < tolv) toler = tolv;
426     }
427
428 //    Edge on same vertex : small one ?
429     if (VA.IsSame(VB)) {
430       Standard_Real cf = 0.,cl = 0.;
431       Handle(Geom_Curve) C3D;
432       if (!BRep_Tool::Degenerated(E)) C3D = BRep_Tool::Curve (E,cf,cl);
433       if (C3D.IsNull()) continue;  // DGNR
434       Standard_Integer np = 0;
435       gp_Pnt deb = C3D->Value(cf);
436       MinMaxPnt (deb,np,minx,miny,minz,maxx,maxy,maxz);
437       gp_Pnt fin = C3D->Value(cl);
438       MinMaxPnt (fin,np,minx,miny,minz,maxx,maxy,maxz);
439       gp_Pnt mid = C3D->Value( (cf+cl)/2. );
440       MinMaxPnt (mid,np,minx,miny,minz,maxx,maxy,maxz);
441       if (!MinMaxSmall (minx,miny,minz,maxx,maxy,maxz,toler)) return Standard_False;
442     } else {
443 //    Other case : two maximum allowed
444       nb ++;
445       if (nb > 2) return Standard_False;
446       if (nb == 1)  {  V1 = VA;  V2 = VB;  E1 = E;  }
447       else if (nb == 2) {
448         if (V1.IsSame(VA) && !V2.IsSame(VB)) return Standard_False;
449         if (V1.IsSame(VB) && !V2.IsSame(VA)) return Standard_False;
450         E2 = E;
451       }
452       else return Standard_False;
453     }
454   }
455
456   if (nb < 2) return Standard_False;   // only one vertex : cannot be a strip ...
457
458 //   Checking if E1 and E2 define a Strip
459   Standard_Real dmax;
460  if (!CheckStripEdges (E1,E2,tol,dmax)) return Standard_False;
461  myStatusStrip = ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
462  return Standard_True;
463 }
464
465 //=======================================================================
466 //function : CheckStripFace
467 //purpose  : 
468 //=======================================================================
469
470 Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckStripFace(const TopoDS_Face& F,
471                                                                TopoDS_Edge& E1, TopoDS_Edge& E2,const Standard_Real tol) 
472 {
473
474   // Standard_Integer stat;
475   if(CheckSingleStrip (F,E1,E2,tol)) return Standard_True ;   // it is a strip
476
477 //    IsStripSupport used as rejection. But this kind of test may be done
478 //    on ANY face, once we are SURE that FindStripEdges is reliable (and fast
479 //    enough)
480
481 //  ?? record a diagnostic StripFace, but without yet lists of edges
482 //  ??  Record Diagnostic "StripFace", no data (should be "Edges1" "Edges2")
483 //      but direction is known (1:U  2:V)
484    // TopoDS_Edge E1,E2;
485     Standard_Real dmax;
486     if(FindStripEdges (F,E1,E2,tol,dmax)) return Standard_True;
487
488 //  Now, trying edges : if there are 2 and only 2 edges greater than tolerance
489 //   (given or sum of vertex tolerances), do they define a strip
490 //  Warning : if yes, they bring different vertices ...
491
492   return Standard_False;
493
494 }
495
496 //=======================================================================
497 //function : CheckSplittingVertices
498 //purpose  : 
499 //=======================================================================
500
501  Standard_Integer ShapeAnalysis_CheckSmallFace::CheckSplittingVertices(const TopoDS_Face& F,
502                                                                        TopTools_DataMapOfShapeListOfShape& MapEdges,
503                                                                        ShapeAnalysis_DataMapOfShapeListOfReal& MapParam,
504                                                                        TopoDS_Compound& theAllVert)
505 {
506
507 //  Prepare array of vertices with their locations //TopTools
508   Standard_Integer nbv = 0, nbp = 0;
509   //TopoDS_Compound theAllVert;
510   BRep_Builder theBuilder;
511   //theBuilder.MakeCompound(theAllVert);
512   TopExp_Explorer itv; // svv Jan11 2000 : porting on DEC
513   for (itv.Init(F,TopAbs_VERTEX); itv.More(); itv.Next()) nbv ++;
514
515   if (nbv == 0) return 0;
516   TopTools_Array1OfShape vtx (1,nbv);
517   TColgp_Array1OfPnt vtp (1,nbv);
518   TColStd_Array1OfReal vto (1,nbv);
519
520   nbp = 0;
521   for (itv.Init(F,TopAbs_VERTEX); itv.More(); itv.Next()) {
522     nbp ++;
523     TopoDS_Vertex unv = TopoDS::Vertex (itv.Current());
524     vtx.SetValue (nbp,unv);
525     gp_Pnt unp = BRep_Tool::Pnt (unv);
526     vtp.SetValue (nbp,unp);
527     Standard_Real unt = myPrecision;
528     if (unt < 0) unt =BRep_Tool::Tolerance (unv);
529     vto.SetValue (nbp,unt);
530   }
531   nbv = nbp;  nbp = 0;  // now, counting splitting vertices
532
533 //  Check edges : are vertices (other than extremities) confused with it ?
534   ShapeAnalysis_Curve SAC;
535   for (Standard_Integer iv = 1; iv <= nbv; iv ++) {
536     TopoDS_Vertex V = TopoDS::Vertex (vtx.Value(iv));
537     TopTools_ListOfShape listEdge;
538     TColStd_ListOfReal listParam;
539     Standard_Boolean issplit = Standard_False;
540     for (TopExp_Explorer ite(F,TopAbs_EDGE); ite.More(); ite.Next()) {
541       TopoDS_Edge E = TopoDS::Edge (ite.Current());
542       TopoDS_Vertex V1,V2;
543       TopExp::Vertices (E,V1,V2);
544       Standard_Real cf,cl;
545       Handle(Geom_Curve) C3D = BRep_Tool::Curve (E,cf,cl);
546       if (C3D.IsNull()) continue;
547       if (V.IsSame(V1) || V.IsSame(V2)) continue;
548       gp_Pnt unp = vtp.Value(iv);
549       Standard_Real unt = vto.Value(iv);
550       gp_Pnt proj;
551       Standard_Real param;
552       Standard_Real dist = SAC.Project (C3D,unp,unt*10.,proj,param,cf,cl);
553       if (dist == 0.0) continue; //smh
554 //  Splitting Vertex to record ?
555       if (dist < unt) {
556 //  If Split occurs at beginning or end, it is not a split ...
557         Standard_Real fpar, lpar, eps = 1.e-06;
558         if (param >=cl || param <= cf) continue; // Out of range
559         fpar = param - cf; lpar = param - cl;
560         if ((Abs(fpar) < eps) || (Abs(lpar) < eps)) continue; // Near end or start 
561         listEdge.Append(E);
562         listParam.Append(param);
563         issplit = Standard_True;
564         
565       }
566     }
567     if(issplit) {
568       nbp ++;
569       theBuilder.Add(theAllVert, V);
570       MapEdges.Bind(V,listEdge);
571       MapParam.Bind(V,listParam);
572     }
573   }
574   if(nbp != 0) 
575   myStatusSplitVert = ShapeExtend::EncodeStatus (ShapeExtend_DONE);
576   return nbp;
577 }
578
579  
580 static Standard_Integer IsoStat
581   (const TColgp_Array2OfPnt& poles,
582    const Standard_Integer uorv, const Standard_Integer rank,
583    const Standard_Real tolpin, const Standard_Real toler)
584 {
585   Standard_Integer i, np = 0;
586   Standard_Integer i0 = (uorv == 1 ? poles.LowerCol() : poles.LowerRow());
587   Standard_Integer i1 = (uorv == 1 ? poles.UpperCol() : poles.UpperRow());
588   Standard_Real xmin = 0.,ymin = 0.,zmin = 0., xmax = 0.,ymax = 0.,zmax = 0.;
589   for (i = i0; i <= i1; i ++) {
590     if (uorv == 1) MinMaxPnt (poles(rank,i),np,xmin,ymin,zmin, xmax,ymax,zmax);
591     else      MinMaxPnt (poles(i,rank), np, xmin,ymin,zmin, xmax,ymax,zmax);
592   }
593   if (MinMaxSmall (xmin,ymin,zmin, xmax,ymax,zmax, tolpin)) return 0;
594   if (MinMaxSmall (xmin,ymin,zmin, xmax,ymax,zmax, toler))  return 1;
595   return 2;
596 }
597
598 static Standard_Boolean CheckPoles(const TColgp_Array2OfPnt& poles, Standard_Integer uorv, Standard_Integer rank)
599 {
600   Standard_Integer i0 = (uorv == 1 ? poles.LowerCol() : poles.LowerRow());
601   Standard_Integer i1 = (uorv == 1 ? poles.UpperCol() : poles.UpperRow());
602   for (Standard_Integer i = i0; i <= i1-1; i ++) {
603     if (uorv == 1) {
604       if(poles(rank,i).IsEqual(poles(rank, i+1), 1e-15)) return Standard_True;
605     } else
606       if(poles(i,rank).IsEqual(poles(i+1,rank), 1e-15)) return Standard_True;
607   }  
608   return Standard_False;
609 }
610 //=======================================================================
611 //function : CheckPin
612 //purpose  : 
613 //=======================================================================
614
615 Standard_Boolean  ShapeAnalysis_CheckSmallFace::CheckPin (const TopoDS_Face& F, Standard_Integer& whatrow,Standard_Integer& sens)
616 {
617   TopLoc_Location loc;
618   Handle(Geom_Surface) surf = BRep_Tool::Surface (F,loc);
619   if (surf->IsKind(STANDARD_TYPE(Geom_ElementarySurface))) return Standard_False;
620
621   Standard_Real toler = myPrecision;
622   if (toler < 0) toler = 1.e-4;
623   Standard_Real tolpin = 1.e-9;  // for sharp sharp pin
624
625 //  Checking the poles
626
627 //  Take the poles : they give good idea of sharpness of a pin
628   Standard_Integer nbu = 0 , nbv = 0;
629   Handle(Geom_BSplineSurface) bs = Handle(Geom_BSplineSurface)::DownCast(surf);
630   Handle(Geom_BezierSurface)  bz = Handle(Geom_BezierSurface)::DownCast(surf);
631   if (!bs.IsNull())  {  nbu = bs->NbUPoles();  nbv = bs->NbVPoles();  }
632   if (!bz.IsNull())  {  nbu = bz->NbUPoles();  nbv = bz->NbVPoles();  }
633   if (nbu == 0 || nbv == 0) return Standard_False;
634
635   TColgp_Array2OfPnt allpoles (1,nbu,1,nbv);
636   if (!bs.IsNull()) bs->Poles (allpoles);
637   if (!bz.IsNull()) bz->Poles (allpoles);
638
639 //  Check each natural bound if it is a singularity (i.e. a pin)
640
641   sens = 0; 
642   Standard_Integer stat = 0;    // 0 none, 1 in U, 2 in V
643   whatrow = 0;  // 0 no row, else rank of row
644   stat = IsoStat(allpoles,1,  1,tolpin,toler);
645   if (stat) { sens = 1;  whatrow = nbu; }
646   
647   stat = IsoStat(allpoles,1,nbu,tolpin,toler);
648   if (stat) { sens = 1; whatrow = nbu;  }
649   
650   stat = IsoStat(allpoles,2,  1,tolpin,toler);
651   if (stat) { sens = 2; whatrow = 1; }
652   
653   stat = IsoStat(allpoles,2,nbv,tolpin,toler);
654   if (stat) { sens = 2; whatrow = nbv;  }
655
656   if (!sens) return Standard_False;    // no pin
657   
658   switch(stat) {
659   case 1: myStatusPin = ShapeExtend::EncodeStatus (ShapeExtend_DONE1); break;
660   case 2: myStatusPin = ShapeExtend::EncodeStatus (ShapeExtend_DONE2); break;
661     default : break;
662   }
663   //  cout<<(whatstat == 1 ? "Smooth" : "Sharp")<<" Pin on "<<(sens == 1 ? "U" : "V")<<" Row n0 "<<whatrow<<endl;
664   if (stat == 1 ) 
665     {
666       // Standard_Boolean EqualPoles = Standard_False;
667       if(CheckPoles(allpoles, 2, nbv)|| CheckPoles(allpoles, 2, 1)
668          ||CheckPoles(allpoles, 1, nbu)|| CheckPoles(allpoles, 1, 1))       
669         myStatusPin = ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
670     }
671   
672   
673   return Standard_True;
674 }
675
676 static Standard_Real TwistedNorm
677 (const Standard_Real x1, const Standard_Real y1, const Standard_Real z1, const Standard_Real x2, const Standard_Real y2, const Standard_Real z2)
678 {  return (x1*x2) + (y1*y2) + (z1*z2);  }
679
680 //=======================================================================
681 //function : CheckTwisted
682 //purpose  : 
683 //=======================================================================
684
685 Standard_Boolean  ShapeAnalysis_CheckSmallFace::CheckTwisted (const TopoDS_Face& F, Standard_Real& paramu,
686                                                              Standard_Real& paramv)
687 {
688   TopLoc_Location loc;
689   Handle(Geom_Surface) surf = BRep_Tool::Surface (F,loc);
690   if (surf->IsKind(STANDARD_TYPE(Geom_ElementarySurface))) return Standard_False;
691
692   Standard_Real toler = myPrecision;
693   if (toler < 0) toler = 1.e-4;
694 ////  GeomLProp_SLProps GLS (surf,2,toler);
695   GeomAdaptor_Surface GAS (surf);
696
697 // to be done : on isos of the surface
698 //  and on edges, at least of outer wire
699   Standard_Integer nbint = 5;
700   TColStd_Array2OfReal nx (1,nbint+1,1,nbint+1);
701   TColStd_Array2OfReal ny (1,nbint+1,1,nbint+1);
702   TColStd_Array2OfReal nz (1,nbint+1,1,nbint+1);
703   Standard_Integer iu,iv;
704   Standard_Real umin,umax,vmin,vmax;
705   surf->Bounds (umin,umax,vmin,vmax);
706   Standard_Real u = umin, du = (umax-umin)/nbint;
707   Standard_Real v = vmin, dv = (umax-umin)/nbint;
708
709   // gp_Dir norm;
710   for (iu = 1; iu <= nbint; iu ++) {
711     for (iv = 1; iv <= nbint; iv ++) {
712 //      GLS.SetParameters (u,v);
713 //      if (GLS.IsNormalDefined()) norm = GLS.Normal();
714       gp_Pnt curp;  gp_Vec V1,V2,VXnorm;
715       GAS.D1 (u,v,curp,V1,V2);
716       VXnorm = V1.Crossed(V2);
717       nx.SetValue (iu,iv,VXnorm.X());
718       ny.SetValue (iu,iv,VXnorm.Y());
719       nz.SetValue (iu,iv,VXnorm.Z());
720       v += dv;
721     }
722     u += du;
723     v = vmin;
724   }
725
726 //  Now, comparing normals on support surface, in both senses
727 //  In principle, it suffuces to check within outer bound
728
729   for (iu = 1; iu < nbint; iu ++) {
730     for (iv = 1; iv < nbint; iv ++) {
731 // We here check each normal (iu,iv) with (iu,iv+1) and with (iu+1,iv)
732 // if for each test, we have negative scalar product, this means angle > 90deg
733 // it is the criterion to say it is twisted
734       if (TwistedNorm ( nx(iu,iv),ny(iu,iv),nz(iu,iv) , nx(iu,iv+1),ny(iu,iv+1),nz(iu,iv+1) ) < 0. ||
735           TwistedNorm ( nx(iu,iv),ny(iu,iv),nz(iu,iv) , nx(iu+1,iv),ny(iu+1,iv),nz(iu+1,iv) ) < 0. ) {
736         myStatusTwisted = ShapeExtend::EncodeStatus (ShapeExtend_DONE);
737         paramu = umin+du*iu-du/2;
738         paramv = vmin+dv*iv-dv/2;
739         return Standard_True;
740       }
741     }
742   }
743
744 //   Now, comparing normals on edges ... to be done
745
746   return Standard_False;
747 }
748
749
750 //=======================================================================
751 //function : CheckPinFace
752 //purpose  : 
753 //=======================================================================
754 // Warning: This function not tested on many examples
755
756  Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckPinFace(const TopoDS_Face& F,
757   TopTools_DataMapOfShapeShape& mapEdges,const Standard_Real toler) 
758 {
759   //ShapeFix_Wire sfw;
760   TopExp_Explorer exp_w (F,TopAbs_WIRE); 
761   exp_w.More();
762   Standard_Real coef1=0, coef2; // =0 for deleting warning (skl)
763   TopoDS_Wire theCurWire = TopoDS::Wire (exp_w.Current());
764   ShapeAnalysis_WireOrder wi;
765   ShapeAnalysis_Wire sfw;
766   Handle(ShapeExtend_WireData) sbwd = new ShapeExtend_WireData(theCurWire);
767   sfw.Load(sbwd);
768   sfw.CheckOrder(wi);
769   Handle(TopTools_HSequenceOfShape) newedges = new TopTools_HSequenceOfShape();
770   Standard_Integer nb = wi.NbEdges();
771   Standard_Integer i = 0;
772   for ( i=1; i <= nb; i++ )
773     newedges->Append ( sbwd->Edge ( wi.Ordered(i) ) );
774   for ( i=1; i <= nb; i++ ) 
775     sbwd->Set ( TopoDS::Edge ( newedges->Value(i) ), i );
776   //sfw.Init(theCurWire,  F, Precision::Confusion());
777   //sfw.FixReorder();
778   //theCurWire = sfw.Wire();
779   theCurWire = sbwd->Wire();
780   i=1;
781   Standard_Boolean done = Standard_False;
782   Standard_Real tol = Precision::Confusion();
783   TopoDS_Edge theFirstEdge, theSecondEdge;
784   Standard_Real d1=0,d2=0;
785   for (TopExp_Explorer exp_e (F,TopAbs_EDGE); exp_e.More(); exp_e.Next()) 
786     {
787       TopoDS_Vertex V1,V2;
788       gp_Pnt p1, p2;
789       if (i==1) 
790         { 
791           theFirstEdge = TopoDS::Edge (exp_e.Current()); 
792           V1 = TopExp::FirstVertex(theFirstEdge);
793           V2 = TopExp::LastVertex(theFirstEdge);
794           p1 = BRep_Tool::Pnt(V1);
795           p2 = BRep_Tool::Pnt(V2);
796           tol = Max(BRep_Tool::Tolerance(V1), BRep_Tool::Tolerance(V2));
797           if (toler > 0) //tol = Max(tol, toler); gka
798           tol = toler;
799           d1 = p1.Distance(p2);
800           if (d1 == 0) return Standard_False;
801           if (d1/tol>=1) coef1 = d1/tol; else continue;
802           if (coef1<=3) continue;
803           i++; 
804           continue;
805         }
806       //Check the length of edge
807       theSecondEdge = TopoDS::Edge (exp_e.Current());
808       V1 = TopExp::FirstVertex(theSecondEdge);
809       V2 = TopExp::LastVertex(theSecondEdge);
810       
811       p1 = BRep_Tool::Pnt(V1);
812       p2 = BRep_Tool::Pnt(V2);
813       if (toler == -1) tol = Max(BRep_Tool::Tolerance(V1), BRep_Tool::Tolerance(V2)); 
814       else tol= toler;
815       if (p1.Distance(p2)> tol) continue;
816       //If there are two pin edges, record them in diagnostic
817       d2 = p1.Distance(p2); //gka
818       if (d2 == 0) return Standard_False;
819       if (d2/tol >= 1) coef2 = d2/tol; else continue;
820       if (coef2<=3) continue;
821       if (coef1>coef2*10) continue;
822       if (coef2>coef1*10)
823         {
824           theFirstEdge = theSecondEdge;
825           coef1 = coef2;
826           continue;
827         }
828       
829       if (CheckPinEdges(theFirstEdge, theSecondEdge, coef1, coef2,toler))
830         {
831           mapEdges.Bind(theFirstEdge,theSecondEdge);
832           myStatusPinFace = ShapeExtend::EncodeStatus (ShapeExtend_DONE);
833           done = Standard_True;
834         }     
835       
836       theFirstEdge = theSecondEdge;
837       coef1 = coef2;
838       //d1 = d2;
839     }
840   return done;
841 }
842
843
844 //=======================================================================
845 //function : CheckPinEdges
846 //purpose  : 
847 //=======================================================================
848 // Warning: This function not tested on many examples
849
850  Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckPinEdges(const TopoDS_Edge& theFirstEdge,const TopoDS_Edge& theSecondEdge,const Standard_Real coef1,
851    const Standard_Real coef2,const Standard_Real toler) const
852 {
853
854  Standard_Real cf1,cl1,cf2,cl2;
855  Handle(Geom_Curve) C1,C2,C3;
856  C1 = BRep_Tool::Curve (theFirstEdge,cf1,cl1);
857  C2 = BRep_Tool::Curve (theSecondEdge,cf2,cl2);
858  gp_Pnt p1, p2, pp1, pp2, pv;
859  Standard_Real d1 = (cf1-cl1)/coef1;
860  Standard_Real d2 = (cf2-cl2)/coef2;
861  //Standard_Real d1 = cf1-cl1/30; //10; gka
862   //Standard_Real d2 = cf2-cl2/30; //10;
863   p1 = C1->Value(cf1);
864   p2 = C1->Value(cl1);
865   pp1 = C2->Value(cf2);
866   pp2 = C2->Value(cl2);
867   Standard_Real tol;
868   Standard_Real paramc1=0, paramc2=0; // =0 for deleting warning (skl)
869   TopoDS_Vertex theSharedV = TopExp::LastVertex(theFirstEdge);
870   if  (toler == -1)  tol = BRep_Tool::Tolerance(theSharedV); else tol = toler;
871   pv = BRep_Tool::Pnt(theSharedV);
872   if (pv.Distance(p1)<=tol) paramc1 = cf1;
873   else if(pv.Distance(p2)<=tol) paramc1 = cl1;
874   if (pv.Distance(pp1)<=tol) paramc2 = cf2;
875   else if(pv.Distance(pp2)<=tol) paramc2 = cl2;
876   //Computing first derivative vectors and compare angle
877 //   gp_Vec V11, V12, V21, V22;
878 //   gp_Pnt tmp;
879 //   C1->D2(paramc1, tmp, V11, V21);
880 //   C2->D2(paramc2, tmp, V12, V22);
881 //   Standard_Real angle1, angle2;
882 //   try{
883 //     angle1 = V11.Angle(V12);
884 //     angle2 = V21.Angle(V22);
885 //   }
886 //   catch (Standard_Failure)
887 //     {
888 //       cout << "Couldn't compute angle between derivative vectors"  <<endl;
889 //       return Standard_False;
890 //     }
891 //   cout << "angle1 "   << angle1<< endl;
892 //   cout << "angle2 "   << angle2<< endl;
893 //   if (angle1<=0.0001) return Standard_True;
894   gp_Pnt proj;
895   if (p1.Distance(p2)<pp1.Distance(pp2)) 
896     {
897       C3=C1;
898       if (paramc1==cf1)
899        proj = C1->Value(paramc1 + (coef1-3)*d1);
900       else proj = C1->Value(paramc1-3*d1);
901         //proj = C1->Value(paramc1 + 9*d1);
902       //else proj = C1->Value(paramc1-d1);
903     }
904   else 
905     {
906       C3=C2;
907       if (paramc2==cf2)
908         proj = C2->Value(paramc2 + (coef2-3)*d2); 
909       else proj = C2->Value(paramc2 -3*d2);
910         //proj = C2->Value(paramc2 + 9*d2); 
911       //else proj = C2->Value(paramc2 -d2); 
912     }
913   Standard_Real param;
914   GeomAdaptor_Curve GAC(C3);
915   Standard_Real f = C3->FirstParameter();
916   Standard_Real l = C3->LastParameter();
917   gp_Pnt result;
918   ShapeAnalysis_Curve SAC;
919   Standard_Real dist = SAC.Project (GAC,proj,tol,result,param);
920   //pdn check if parameter of projection is in the domain of the edge.
921   if (param < f || param > l) return Standard_False;
922   if (dist > tol) return Standard_False;
923   if (dist <= tol) {
924      //Computing first derivative vectors and compare angle
925       gp_Vec V11, V12, V21, V22;
926       gp_Pnt tmp;
927       C1->D2(paramc1, tmp, V11, V21);
928       C2->D2(paramc2, tmp, V12, V22);
929       Standard_Real angle1=0, angle2=0;
930       try{
931         angle1 = V11.Angle(V12);
932         angle2 = V21.Angle(V22);
933       }
934       catch (Standard_Failure)
935         {
936           cout << "Couldn't compute angle between derivative vectors"  <<endl;
937           return Standard_False;
938         }
939 //       cout << "angle1 "   << angle1<< endl;
940 //       cout << "angle2 "   << angle2<< endl;
941       if ((angle1<=0.001 && angle2<=0.01) || ((M_PI-angle2)<= 0.001 && (M_PI-angle2)<= 0.01)) return Standard_True;
942       else return Standard_False;
943     } 
944     
945   return Standard_False;
946 }
947