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) if(poles(rank,i).IsEqual(poles(rank, i+1), 1e-15)) return Standard_True;
 586     else      if(poles(i,rank).IsEqual(poles(i+1,rank), 1e-15)) return Standard_True;
 587   }
 588   return Standard_False;
 589 }
 590 //=======================================================================
 591 //function : CheckPin
 592 //purpose  :
 593 //=======================================================================
 594
 595 Standard_Boolean  ShapeAnalysis_CheckSmallFace::CheckPin (const TopoDS_Face& F, Standard_Integer& whatrow,Standard_Integer& sens)
 596 {
 597   TopLoc_Location loc;
 598   Handle(Geom_Surface) surf = BRep_Tool::Surface (F,loc);
 599   if (surf->IsKind(STANDARD_TYPE(Geom_ElementarySurface))) return Standard_False;
 600
 601   Standard_Real toler = myPrecision;
 602   if (toler < 0) toler = 1.e-4;
 603   Standard_Real tolpin = 1.e-9;  // for sharp sharp pin
 604
 605 //  Checking the poles
 606
 607 //  Take the poles : they give good idea of sharpness of a pin
 608   Standard_Integer nbu = 0 , nbv = 0;
 609   Handle(Geom_BSplineSurface) bs = Handle(Geom_BSplineSurface)::DownCast(surf);
 610   Handle(Geom_BezierSurface)  bz = Handle(Geom_BezierSurface)::DownCast(surf);
 611   if (!bs.IsNull())  {  nbu = bs->NbUPoles();  nbv = bs->NbVPoles();  }
 612   if (!bz.IsNull())  {  nbu = bz->NbUPoles();  nbv = bz->NbVPoles();  }
 613   if (nbu == 0 || nbv == 0) return Standard_False;
 614
 615   TColgp_Array2OfPnt allpoles (1,nbu,1,nbv);
 616   if (!bs.IsNull()) bs->Poles (allpoles);
 617   if (!bz.IsNull()) bz->Poles (allpoles);
 618
 619 //  Check each natural bound if it is a singularity (i.e. a pin)
 620
 621   sens = 0;
 622   Standard_Integer stat = 0;    // 0 none, 1 in U, 2 in V
 623   whatrow = 0;  // 0 no row, else rank of row
 624   stat = IsoStat(allpoles,1,  1,tolpin,toler);
 625   if (stat) { sens = 1;  whatrow = nbu; }
 626
 627   stat = IsoStat(allpoles,1,nbu,tolpin,toler);
 628   if (stat) { sens = 1; whatrow = nbu;  }
 629
 630   stat = IsoStat(allpoles,2,  1,tolpin,toler);
 631   if (stat) { sens = 2; whatrow = 1; }
 632
 633   stat = IsoStat(allpoles,2,nbv,tolpin,toler);
 634   if (stat) { sens = 2; whatrow = nbv;  }
 635
 636   if (!sens) return Standard_False;    // no pin
 637
 638   switch(stat) {
 639   case 1: myStatusPin = ShapeExtend::EncodeStatus (ShapeExtend_DONE1); break;
 640   case 2: myStatusPin = ShapeExtend::EncodeStatus (ShapeExtend_DONE2); break;
 641     default : break;
 642   }
 643   //  cout<<(whatstat == 1 ? "Smooth" : "Sharp")<<" Pin on "<<(sens == 1 ? "U" : "V")<<" Row n0 "<<whatrow<<endl;
 644   if (stat == 1 )
 645     {
 646       // Standard_Boolean EqualPoles = Standard_False;
 647       if(CheckPoles(allpoles, 2, nbv)|| CheckPoles(allpoles, 2, 1)
 648          ||CheckPoles(allpoles, 1, nbu)|| CheckPoles(allpoles, 1, 1))
 649         myStatusPin = ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
 650     }
 651
 652
 653   return Standard_True;
 654 }
 655
 656 static Standard_Real TwistedNorm
 657 (const Standard_Real x1, const Standard_Real y1, const Standard_Real z1, const Standard_Real x2, const Standard_Real y2, const Standard_Real z2)
 658 {  return (x1*x2) + (y1*y2) + (z1*z2);  }
 659
 660 //=======================================================================
 661 //function : CheckTwisted
 662 //purpose  :
 663 //=======================================================================
 664
 665 Standard_Boolean  ShapeAnalysis_CheckSmallFace::CheckTwisted (const TopoDS_Face& F, Standard_Real& paramu,
 666                                                              Standard_Real& paramv)
 667 {
 668   TopLoc_Location loc;
 669   Handle(Geom_Surface) surf = BRep_Tool::Surface (F,loc);
 670   if (surf->IsKind(STANDARD_TYPE(Geom_ElementarySurface))) return Standard_False;
 671
 672   Standard_Real toler = myPrecision;
 673   if (toler < 0) toler = 1.e-4;
 674 ////  GeomLProp_SLProps GLS (surf,2,toler);
 675   GeomAdaptor_Surface GAS (surf);
 676
 677 // to be done : on isos of the surface
 678 //  and on edges, at least of outer wire
 679   Standard_Integer nbint = 5;
 680   TColStd_Array2OfReal nx (1,nbint+1,1,nbint+1);
 681   TColStd_Array2OfReal ny (1,nbint+1,1,nbint+1);
 682   TColStd_Array2OfReal nz (1,nbint+1,1,nbint+1);
 683   Standard_Integer iu,iv;
 684   Standard_Real umin,umax,vmin,vmax;
 685   surf->Bounds (umin,umax,vmin,vmax);
 686   Standard_Real u = umin, du = (umax-umin)/nbint;
 687   Standard_Real v = vmin, dv = (umax-umin)/nbint;
 688
 689   // gp_Dir norm;
 690   for (iu = 1; iu <= nbint; iu ++) {
 691     for (iv = 1; iv <= nbint; iv ++) {
 692 //      GLS.SetParameters (u,v);
 693 //      if (GLS.IsNormalDefined()) norm = GLS.Normal();
 694       gp_Pnt curp;  gp_Vec V1,V2,VXnorm;
 695       GAS.D1 (u,v,curp,V1,V2);
 696       VXnorm = V1.Crossed(V2);
 697       nx.SetValue (iu,iv,VXnorm.X());
 698       ny.SetValue (iu,iv,VXnorm.Y());
 699       nz.SetValue (iu,iv,VXnorm.Z());
 700       v += dv;
 701     }
 702     u += du;
 703     v = vmin;
 704   }
 705
 706 //  Now, comparing normals on support surface, in both senses
 707 //  In principle, it suffuces to check within outer bound
 708
 709   for (iu = 1; iu < nbint; iu ++) {
 710     for (iv = 1; iv < nbint; iv ++) {
 711 // We here check each normal (iu,iv) with (iu,iv+1) and with (iu+1,iv)
 712 // if for each test, we have negative scalar product, this means angle > 90deg
 713 // it is the criterion to say it is twisted
 714       if (TwistedNorm ( nx(iu,iv),ny(iu,iv),nz(iu,iv) , nx(iu,iv+1),ny(iu,iv+1),nz(iu,iv+1) ) < 0. ||
 715           TwistedNorm ( nx(iu,iv),ny(iu,iv),nz(iu,iv) , nx(iu+1,iv),ny(iu+1,iv),nz(iu+1,iv) ) < 0. ) {
 716         myStatusTwisted = ShapeExtend::EncodeStatus (ShapeExtend_DONE);
 717         paramu = umin+du*iu-du/2;
 718         paramv = vmin+dv*iv-dv/2;
 719         return Standard_True;
 720       }
 721     }
 722   }
 723
 724 //   Now, comparing normals on edges ... to be done
 725
 726   return Standard_False;
 727 }
 728
 729
 730 //=======================================================================
 731 //function : CheckPinFace
 732 //purpose  :
 733 //=======================================================================
 734 // Warning: This function not tested on many examples
 735
 736  Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckPinFace(const TopoDS_Face& F,
 737   TopTools_DataMapOfShapeShape& mapEdges,const Standard_Real toler)
 738 {
 739   //ShapeFix_Wire sfw;
 740   TopExp_Explorer exp_w (F,TopAbs_WIRE);
 741   exp_w.More();
 742   Standard_Real coef1=0, coef2; // =0 for deleting warning (skl)
 743   TopoDS_Wire theCurWire = TopoDS::Wire (exp_w.Current());
 744   ShapeAnalysis_WireOrder wi;
 745   ShapeAnalysis_Wire sfw;
 746   Handle(ShapeExtend_WireData) sbwd = new ShapeExtend_WireData(theCurWire);
 747   sfw.Load(sbwd);
 748   sfw.CheckOrder(wi);
 749   Handle(TopTools_HSequenceOfShape) newedges = new TopTools_HSequenceOfShape();
 750   Standard_Integer nb = wi.NbEdges();
 751   Standard_Integer i = 0;
 752   for ( i=1; i <= nb; i++ )
 753     newedges->Append ( sbwd->Edge ( wi.Ordered(i) ) );
 754   for ( i=1; i <= nb; i++ )
 755     sbwd->Set ( TopoDS::Edge ( newedges->Value(i) ), i );
 756   //sfw.Init(theCurWire,  F, Precision::Confusion());
 757   //sfw.FixReorder();
 758   //theCurWire = sfw.Wire();
 759   theCurWire = sbwd->Wire();
 760   i=1;
 761   Standard_Boolean done = Standard_False;
 762   Standard_Real tol = Precision::Confusion();
 763   TopoDS_Edge theFirstEdge, theSecondEdge;
 764   Standard_Real d1=0,d2=0;
 765   for (TopExp_Explorer exp_e (F,TopAbs_EDGE); exp_e.More(); exp_e.Next())
 766     {
 767       TopoDS_Vertex V1,V2;
 768       gp_Pnt p1, p2;
 769       if (i==1)
 770         {
 771           theFirstEdge = TopoDS::Edge (exp_e.Current());
 772           V1 = TopExp::FirstVertex(theFirstEdge);
 773           V2 = TopExp::LastVertex(theFirstEdge);
 774           p1 = BRep_Tool::Pnt(V1);
 775           p2 = BRep_Tool::Pnt(V2);
 776           tol = Max(BRep_Tool::Tolerance(V1), BRep_Tool::Tolerance(V2));
 777           if (toler > 0) //tol = Max(tol, toler); gka
 778           tol = toler;
 779           d1 = p1.Distance(p2);
 780           if (d1 == 0) return Standard_False;
 781           if (d1/tol>=1) coef1 = d1/tol; else continue;
 782           if (coef1<=3) continue;
 783           i++;
 784           continue;
 785         }
 786       //Check the length of edge
 787       theSecondEdge = TopoDS::Edge (exp_e.Current());
 788       V1 = TopExp::FirstVertex(theSecondEdge);
 789       V2 = TopExp::LastVertex(theSecondEdge);
 790
 791       p1 = BRep_Tool::Pnt(V1);
 792       p2 = BRep_Tool::Pnt(V2);
 793       if (toler == -1) tol = Max(BRep_Tool::Tolerance(V1), BRep_Tool::Tolerance(V2));
 794       else tol= toler;
 795       if (p1.Distance(p2)> tol) continue;
 796       //If there are two pin edges, record them in diagnostic
 797       d2 = p1.Distance(p2); //gka
 798       if (d2 == 0) return Standard_False;
 799       if (d2/tol >= 1) coef2 = d2/tol; else continue;
 800       if (coef2<=3) continue;
 801       if (coef1>coef2*10) continue;
 802       if (coef2>coef1*10)
 803         {
 804           theFirstEdge = theSecondEdge;
 805           coef1 = coef2;
 806           continue;
 807         }
 808
 809       if (CheckPinEdges(theFirstEdge, theSecondEdge, coef1, coef2,toler))
 810         {
 811           mapEdges.Bind(theFirstEdge,theSecondEdge);
 812           myStatusPinFace = ShapeExtend::EncodeStatus (ShapeExtend_DONE);
 813           done = Standard_True;
 814         }
 815
 816       theFirstEdge = theSecondEdge;
 817       coef1 = coef2;
 818       //d1 = d2;
 819     }
 820   return done;
 821 }
 822
 823
 824 //=======================================================================
 825 //function : CheckPinEdges
 826 //purpose  :
 827 //=======================================================================
 828 // Warning: This function not tested on many examples
 829
 830  Standard_Boolean ShapeAnalysis_CheckSmallFace::CheckPinEdges(const TopoDS_Edge& theFirstEdge,const TopoDS_Edge& theSecondEdge,const Standard_Real coef1,
 831    const Standard_Real coef2,const Standard_Real toler) const
 832 {
 833
 834  Standard_Real cf1,cl1,cf2,cl2;
 835  Handle(Geom_Curve) C1,C2,C3;
 836  C1 = BRep_Tool::Curve (theFirstEdge,cf1,cl1);
 837  C2 = BRep_Tool::Curve (theSecondEdge,cf2,cl2);
 838  gp_Pnt p1, p2, pp1, pp2, pv;
 839  Standard_Real d1 = (cf1-cl1)/coef1;
 840  Standard_Real d2 = (cf2-cl2)/coef2;
 841  //Standard_Real d1 = cf1-cl1/30; //10; gka
 842   //Standard_Real d2 = cf2-cl2/30; //10;
 843   p1 = C1->Value(cf1);
 844   p2 = C1->Value(cl1);
 845   pp1 = C2->Value(cf2);
 846   pp2 = C2->Value(cl2);
 847   Standard_Real tol;
 848   Standard_Real paramc1=0, paramc2=0; // =0 for deleting warning (skl)
 849   TopoDS_Vertex theSharedV = TopExp::LastVertex(theFirstEdge);
 850   if  (toler == -1)  tol = BRep_Tool::Tolerance(theSharedV); else tol = toler;
 851   pv = BRep_Tool::Pnt(theSharedV);
 852   if (pv.Distance(p1)<=tol) paramc1 = cf1;
 853   else if(pv.Distance(p2)<=tol) paramc1 = cl1;
 854   if (pv.Distance(pp1)<=tol) paramc2 = cf2;
 855   else if(pv.Distance(pp2)<=tol) paramc2 = cl2;
 856   //Computing first derivative vectors and compare angle
 857 //   gp_Vec V11, V12, V21, V22;
 858 //   gp_Pnt tmp;
 859 //   C1->D2(paramc1, tmp, V11, V21);
 860 //   C2->D2(paramc2, tmp, V12, V22);
 861 //   Standard_Real angle1, angle2;
 862 //   try{
 863 //     angle1 = V11.Angle(V12);
 864 //     angle2 = V21.Angle(V22);
 865 //   }
 866 //   catch (Standard_Failure)
 867 //     {
 868 //       cout << "Couldn't compute angle between derivative vectors"  <<endl;
 869 //       return Standard_False;
 870 //     }
 871 //   cout << "angle1 "   << angle1<< endl;
 872 //   cout << "angle2 "   << angle2<< endl;
 873 //   if (angle1<=0.0001) return Standard_True;
 874   gp_Pnt proj;
 875   if (p1.Distance(p2)<pp1.Distance(pp2))
 876     {
 877       C3=C1;
 878       if (paramc1==cf1)
 879        proj = C1->Value(paramc1 + (coef1-3)*d1);
 880       else proj = C1->Value(paramc1-3*d1);
 881         //proj = C1->Value(paramc1 + 9*d1);
 882       //else proj = C1->Value(paramc1-d1);
 883     }
 884   else
 885     {
 886       C3=C2;
 887       if (paramc2==cf2)
 888         proj = C2->Value(paramc2 + (coef2-3)*d2);
 889       else proj = C2->Value(paramc2 -3*d2);
 890         //proj = C2->Value(paramc2 + 9*d2);
 891       //else proj = C2->Value(paramc2 -d2);
 892     }
 893   Standard_Real param;
 894   GeomAdaptor_Curve GAC(C3);
 895   Standard_Real f = C3->FirstParameter();
 896   Standard_Real l = C3->LastParameter();
 897   gp_Pnt result;
 898   ShapeAnalysis_Curve SAC;
 899   Standard_Real dist = SAC.Project (GAC,proj,tol,result,param);
 900   //pdn check if parameter of projection is in the domain of the edge.
 901   if (param < f || param > l) return Standard_False;
 902   if (dist > tol) return Standard_False;
 903   if (dist <= tol) {
 904      //Computing first derivative vectors and compare angle
 905       gp_Vec V11, V12, V21, V22;
 906       gp_Pnt tmp;
 907       C1->D2(paramc1, tmp, V11, V21);
 908       C2->D2(paramc2, tmp, V12, V22);
 909       Standard_Real angle1=0, angle2=0;
 910       try{
 911         angle1 = V11.Angle(V12);
 912         angle2 = V21.Angle(V22);
 913       }
 914       catch (Standard_Failure)
 915         {
 916           cout << "Couldn't compute angle between derivative vectors"  <<endl;
 917           return Standard_False;
 918         }
 919 //       cout << "angle1 "   << angle1<< endl;
 920 //       cout << "angle2 "   << angle2<< endl;
 921       if ((angle1<=0.001 && angle2<=0.01) || ((M_PI-angle2)<= 0.001 && (M_PI-angle2)<= 0.01)) return Standard_True;
 922       else return Standard_False;
 923     }
 924
 925   return Standard_False;
 926 }
 927