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