src/BRepMesh/BRepMesh_Classifier.cxx

   1 // File:        BRepMesh_Classifier.cxx
   2 // Created:     Thu Jun 26 14:54:02 1997
   3 // Author:      Laurent PAINNOT
   4 //              <lpa@penox.paris1.matra-dtv.fr>
   5
   6
   7 #include <BRepMesh_Classifier.ixx>
   8
   9 // Kernel
  10 #include <Precision.hxx>
  11 #include <Standard_ErrorHandler.hxx>
  12 #include <TColStd_ListOfTransient.hxx>
  13 #include <TColStd_Array1OfInteger.hxx>
  14 #include <TColStd_DataMapOfIntegerInteger.hxx>
  15 #include <ElCLib.hxx>
  16 // Geometry
  17 #include <gp_Pnt.hxx>
  18 #include <gp_Pnt2d.hxx>
  19 #include <TColgp_SequenceOfPnt2d.hxx>
  20 #include <TColgp_Array1OfPnt2d.hxx>
  21 #include <GeomAbs_SurfaceType.hxx>
  22 #include <Geom2dInt_Geom2dCurveTool.hxx>
  23 #include <Geom2d_Line.hxx>
  24 #include <Geom2d_BezierCurve.hxx>
  25 #include <Geom2d_BSplineCurve.hxx>
  26 #include <Geom2d_TrimmedCurve.hxx>
  27 // Topology
  28 #include <BRep_Tool.hxx>
  29 #include <BRepTools.hxx>
  30 #include <BRepTools_WireExplorer.hxx>
  31 #include <BRepAdaptor_Curve2d.hxx>
  32 #include <TopAbs_Orientation.hxx>
  33 #include <TopExp.hxx>
  34 #include <TopExp_Explorer.hxx>
  35 #include <TopoDS_Edge.hxx>
  36 #include <TopoDS.hxx>
  37 #include <CSLib_Class2d.hxx>
  38 #include <Poly_PolygonOnTriangulation.hxx>
  39 // BRepMesh
  40 #include <BRepMesh_Vertex.hxx>
  41 #include <BRepMesh_Array1OfBiPoint.hxx>
  42 #include <BRepMesh_PairOfPolygon.hxx>
  43
  44 #ifdef DEB_MESH
  45 static Standard_Integer debwire;
  46 static Standard_Integer debedge;
  47 static Standard_Integer debclass = 0;
  48 #endif
  49 static const Standard_Real MIN_DIST = 2.E-5; //EPA: real mesh is created in the grid 10E5x10E5, so intersection should be cheched
  50                                              //     with double of discretization.
  51 static const Standard_Real PARALL_COND = Sin(PI/3.0);
  52 static const Standard_Real RESOLUTION = 1.0E-16; //OCC319
  53
  54
  55 //=======================================================================
  56 //function : IsLine
  57 //purpose  :
  58 //=======================================================================
  59
  60 static Standard_Boolean IsLine(const Handle(Geom2d_Curve)& C2d)
  61 {
  62   Standard_Boolean IsALine = Standard_False;
  63   if ( C2d->IsKind(STANDARD_TYPE(Geom2d_Line)) )
  64   {
  65     IsALine = Standard_True;
  66   }
  67   else if ( C2d->IsKind(STANDARD_TYPE(Geom2d_BSplineCurve)) )
  68   {
  69     Handle(Geom2d_BSplineCurve) BS = *((Handle(Geom2d_BSplineCurve)*)&C2d);
  70     IsALine = (BS->NbPoles() == 2);
  71   }
  72   else if ( C2d->IsKind(STANDARD_TYPE(Geom2d_BezierCurve)) )
  73   {
  74     Handle(Geom2d_BezierCurve) Bz = *((Handle(Geom2d_BezierCurve)*)&C2d);
  75     IsALine = (Bz->NbPoles() == 2);
  76   }
  77   else if ( C2d->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)) )
  78   {
  79     Handle(Geom2d_TrimmedCurve) Curv = *((Handle(Geom2d_TrimmedCurve)*)&C2d);
  80     IsALine = IsLine(Curv->BasisCurve());
  81   }
  82   return IsALine;
  83 }
  84
  85 //=======================================================================
  86 //function : AnalizeWire
  87 //purpose  :
  88 //=======================================================================
  89
  90 void BRepMesh_Classifier::AnalizeWire (const TColgp_SequenceOfPnt2d&  theSeqPnt2d,
  91                                        const Standard_Real Umin,  const Standard_Real Umax,
  92                                        const Standard_Real Vmin,  const Standard_Real Vmax)
  93 {
  94   const Standard_Integer nbpnts = theSeqPnt2d.Length();
  95   if (nbpnts < 2) return;
  96
  97   // Accumulate angle
  98   TColgp_Array1OfPnt2d PClass(1,nbpnts);
  99   Standard_Integer i, ii;
 100   Standard_Real theangle = 0.0;
 101   gp_Pnt2d p1 = theSeqPnt2d(1), p2 = theSeqPnt2d(2), p3;
 102   PClass(1) = p1;
 103   PClass(2) = p2;
 104   for (i = 1; i <= nbpnts; i++)
 105   {
 106     ii = i + 2;
 107     if (ii > nbpnts)
 108     {
 109       p3 = PClass(ii-nbpnts);
 110     }
 111     else
 112     {
 113       p3 = theSeqPnt2d.Value(ii);
 114       PClass(ii) = p3;
 115     }
 116     gp_Vec2d A(p1,p2), B(p2,p3);
 117     if (A.SquareMagnitude() > 1.e-16 && B.SquareMagnitude() > 1.e-16)
 118     {
 119       const Standard_Real a = A.Angle(B);
 120       const Standard_Real aa = Abs(a);
 121       // Check if vectors are opposite
 122       if (aa > Precision::Angular() && (PI - aa) > Precision::Angular())
 123       {
 124         theangle += a;
 125         p1 = p2;
 126       }
 127     }
 128     p2 = p3;
 129   }
 130   // Check for zero angle - treat self intersecting wire as outer
 131   if (Abs(theangle) < Precision::Angular()) theangle = 0.0;
 132
 133   TabClass.Append((void *)new CSLib_Class2d(PClass,Toluv,Toluv,Umin,Vmin,Umax,Vmax));
 134   TabOrien.Append((theangle < 0.0) ? 0 : 1);
 135 }
 136
 137 //=======================================================================
 138 //function : BRepMesh_Classifier
 139 //purpose  :
 140 //=======================================================================
 141
 142 BRepMesh_Classifier::BRepMesh_Classifier(const TopoDS_Face& aFace,
 143                                          const Standard_Real TolUV,
 144                                          const BRepMesh_DataMapOfShapePairOfPolygon& edges,
 145                                          const TColStd_IndexedMapOfInteger& themap,
 146                                          const Handle(BRepMesh_DataStructureOfDelaun)& Str,
 147                                          const Standard_Real Umin,
 148                                          const Standard_Real Umax,
 149                                          const Standard_Real Vmin,
 150                                          const Standard_Real Vmax):
 151                                          Toluv(TolUV), Face(aFace),
 152                                          myState(BRepMesh_NoError)
 153 {
 154   //-- impasse sur les surfs definies sur plus d une periode
 155
 156   //-- once definition
 157   Face.Orientation(TopAbs_FORWARD);
 158
 159   TopoDS_Edge  edge;
 160   BRepTools_WireExplorer WireExplorer;
 161   //TopExp_Explorer FaceExplorer;
 162   TopoDS_Iterator FaceExplorer;
 163
 164   TColgp_SequenceOfPnt2d aWirePoints, aWire;
 165   TColStd_SequenceOfInteger aWireLength;
 166
 167
 168   //-- twice definitions
 169   TopAbs_Orientation anOr = TopAbs_FORWARD;
 170   Standard_Boolean falsewire = Standard_False;
 171   Standard_Integer i, index, firstindex = 0, lastindex = 0, nbedges = 0;
 172 #ifdef DEB_MESH
 173   debwire = 0;
 174 #endif
 175
 176   for(FaceExplorer.Initialize(Face); FaceExplorer.More(); FaceExplorer.Next())
 177   {
 178     if( FaceExplorer.Value().ShapeType()!= TopAbs_WIRE)
 179       continue;
 180 #ifdef DEB_MESH
 181     if (debclass) { debwire++;  cout <<endl;  cout << "#wire no "<<debwire; debedge = 0;}
 182 #endif
 183     // For each wire we create a data map, linking vertices (only
 184     // the ends of edges) with their positions in the sequence of
 185     // all 2d points from this wire.
 186     // When we meet some vertex for the second time - the piece
 187     // of sequence is treated for a HOLE and quits the sequence.
 188     // Actually, we must unbind the vertices belonging to the
 189     // loop from the map, but since they can't appear twice on the
 190     // valid wire, leave them for a little speed up.
 191     nbedges = 0;
 192     TColgp_SequenceOfPnt2d SeqPnt2d;
 193     TColStd_DataMapOfIntegerInteger NodeInSeq;
 194     // Start traversing the wire
 195     for (WireExplorer.Init(TopoDS::Wire(FaceExplorer.Value()),Face); WireExplorer.More(); WireExplorer.Next())
 196     {
 197       edge = WireExplorer.Current();
 198 #ifdef DEB_MESH
 199       if (debclass) { debedge++; cout << endl; cout << "#edge no "<<debedge <<endl;}
 200 #endif
 201       anOr = edge.Orientation();
 202       if (anOr != TopAbs_FORWARD && anOr != TopAbs_REVERSED) continue;
 203       if (edges.IsBound(edge))
 204       {
 205         // Retrieve polygon
 206         // Define the direction for adding points to SeqPnt2d
 207         Standard_Integer iFirst,iLast,iIncr;
 208         const BRepMesh_PairOfPolygon& pair = edges.Find(edge);
 209         Handle(Poly_PolygonOnTriangulation) NOD;
 210         if (anOr == TopAbs_FORWARD)
 211         {
 212           NOD = pair.First();
 213           iFirst = 1;
 214           iLast  = NOD->NbNodes();
 215           iIncr  = 1;
 216         }
 217         else
 218         {
 219           NOD = pair.Last();
 220           iFirst = NOD->NbNodes();
 221           iLast  = 1;
 222           iIncr  = -1;
 223         }
 224         const TColStd_Array1OfInteger& indices = NOD->Nodes();
 225
 226         // indexFirst and nodeLast are the indices of first and last
 227         // vertices of the edge in IndexedMap <Str>
 228         const Standard_Integer indexFirst = themap.FindKey(indices(iFirst));
 229         const Standard_Integer indexLast = themap.FindKey(indices(iLast));
 230
 231         // Skip degenerated edge : OCC481(apo)
 232         if (indexLast == indexFirst && (iLast-iFirst) == iIncr) continue;
 233
 234         // If there's a gap between edges -> raise <falsewire> flag
 235         if (nbedges)
 236         {
 237           if (indexFirst != lastindex)
 238           {
 239             falsewire = Standard_True;
 240             break;
 241           }
 242         }
 243         else firstindex = indexFirst;
 244               lastindex = indexLast;
 245
 246         // Record first vertex (to detect loops)
 247         NodeInSeq.Bind(indexFirst,SeqPnt2d.Length()+1);
 248
 249         // Add vertices in sequence
 250         for (i = iFirst; i != iLast; i += iIncr)
 251         {
 252           index = (i == iFirst)? indexFirst : themap.FindKey(indices(i));
 253
 254           gp_Pnt2d vp(Str->GetNode(index).Coord());
 255           SeqPnt2d.Append(vp);
 256 #ifdef DEB_MESH
 257           if (debclass) cout<<"point p"<<index<<" "<<vp.X()<<" "<< vp.Y()<<endl;
 258 #endif
 259         }
 260
 261         // Now, is there a loop?
 262         if (NodeInSeq.IsBound(indexLast))
 263         {
 264           // Yes, treat it separately as a hole
 265           // 1. Divide points into main wire and a loop
 266           const Standard_Integer iWireStart = NodeInSeq(indexLast);
 267           if(iWireStart < SeqPnt2d.Length()) {
 268             SeqPnt2d.Split(iWireStart, aWire);
 269             //OCC319->  the operation will be done later
 270             // 2. Proceed the loop
 271             //AnalizeWire(aLoop, Umin, Umax, Vmin, Vmax, aWirePoints, aWireLength, NbBiPoint);
 272             aWireLength.Append(aWire.Length());
 273             aWirePoints.Append(aWire);
 274             //<-OCC319
 275           }
 276         }
 277         nbedges++;
 278       }
 279     }
 280
 281     if (nbedges)
 282     {
 283       // Isn't it open?
 284       if (falsewire || (firstindex != lastindex) || SeqPnt2d.Length() > 1)
 285       {
 286         myState = BRepMesh_OpenWire;
 287         return;
 288       }
 289     }
 290     else
 291     {
 292 #ifdef DEB_MESH
 293       cout <<"Warning : empty wire" <<endl;
 294 #endif
 295     }
 296   }
 297
 298   const Standard_Integer nbwires = aWireLength.Length();
 299   Standard_Integer NbBiPoint = aWirePoints.Length();
 300   BRepMesh_Array1OfBiPoint BiPoints(0,NbBiPoint);
 301
 302   BRepMesh_BiPoint *BP;
 303   Standard_Real *Coordinates1;
 304   Standard_Real x1, y1, x2, y2, xstart, ystart;
 305   Standard_Integer j, l = 1;
 306   BP = &(BiPoints.ChangeValue(1));
 307
 308   // Fill array of segments (bi-points)
 309   for (i = 1; i <= nbwires; i++)
 310   {
 311     const Standard_Integer len = aWireLength(i) + 1;
 312     for (j = 1; j <= len; j++)
 313     {
 314       // Obtain last point of the segment
 315       if (j == len)
 316       {
 317         x2 = xstart;
 318         y2 = ystart;
 319       }
 320       else
 321       {
 322         const gp_Pnt2d& PT = aWirePoints(l); l++;
 323         x2 = PT.X();
 324         y2 = PT.Y();
 325       }
 326       // Build segment (bi-point)
 327       if (j == 1)
 328       {
 329         xstart = x2;
 330         ystart = y2;
 331       }
 332       else
 333       {
 334         Coordinates1 = ((Standard_Real*)(BP->Coordinates())); BP++;
 335         Coordinates1[0] = x1;
 336         Coordinates1[1] = y1;
 337         Coordinates1[2] = x2;
 338         Coordinates1[3] = y2;
 339         Coordinates1[4] = x2 - x1;
 340         Coordinates1[5] = y2 - y1;
 341       }
 342       x1 = x2;
 343       y1 = y2;
 344     }
 345   }
 346
 347   Standard_Real *Coordinates2;
 348   Standard_Real A1, B1, C1, A2, B2, C2, AB, BC, CA, xc, yc;
 349   Standard_Real  mu1, d, mu2;
 350   Standard_Integer ik, ikEnd = 0, jk, jkEnd;
 351   Standard_Real x11, x12, y11, y12, x21, x22, y21, y22;
 352   for(i = 1; i <= nbwires; i++)
 353   {
 354     ik = ikEnd + 1;  ikEnd += aWireLength(i);
 355     // Explore first wire
 356     for (; ik <= ikEnd; ik++)
 357     {
 358       Coordinates1 = ((Standard_Real*)(BiPoints.ChangeValue(ik).Coordinates()));
 359       x11 = Coordinates1[0];
 360       y11 = Coordinates1[1];
 361       x12 = Coordinates1[2];
 362       y12 = Coordinates1[3];
 363       A1 =  Coordinates1[5];
 364       B1 = -Coordinates1[4];
 365       C1 = - x11*A1 - y11*B1;
 366       //mu1 = Sqrt(A1*A1+B1*B1);
 367       mu1 = A1*A1+B1*B1;
 368       for (j = i; j <= nbwires; j++)
 369       {
 370         //for i==j the algorithm check current wire on selfintersection
 371         if (j == i)
 372         {
 373           jk = ik + 2;  jkEnd = ikEnd;
 374         }
 375         else
 376         {
 377           jk = jkEnd + 1;  jkEnd = jk + aWireLength(j) - 1;
 378         }
 379         // Explore second wire
 380         for (; jk <= jkEnd; jk++)
 381         {
 382           // don't check end's segment of the wire on selfrestriction
 383           if (jk == ikEnd) continue;
 384           Coordinates2 = ((Standard_Real*)(BiPoints.ChangeValue(jk).Coordinates()));
 385           x21 = Coordinates2[0];
 386           y21 = Coordinates2[1];
 387           x22 = Coordinates2[2];
 388           y22 = Coordinates2[3];
 389           A2 =  Coordinates2[5];
 390           B2 = -Coordinates2[4];
 391           C2 = - x21*A2 - y21*B2;
 392           //mu2 = Sqrt(A2*A2+B2*B2);
 393           mu2 = A2*A2+B2*B2;
 394           //different segments may have common vertex (see OCC287 bug for example)
 395           //if(x22 == x11 && y22 == y11){ myState = BRepMesh_OpenWire;  return;}
 396           AB = A1*B2 - A2*B1;
 397           //check on minimal of distance between current segment and points of another linear segments - OCC319
 398           //d = Abs(A1*x22 + B1*y22 + C1);
 399           d = A1*x22 + B1*y22 + C1;
 400           if(i != j &&                        // if compared wires are different &&
 401              AB*AB > PARALL_COND*PARALL_COND*mu1*mu2 && // angle between two segments greater then PARALL_COND &&
 402              d*d < MIN_DIST*MIN_DIST*mu1 &&              // distance between vertex of the segment and other one's less then MIN_DIST
 403              (x22-x11)*(x22-x12) < 0.0 && (y22-y11)*(y22-y12) < 0.0)
 404           {
 405             myState = BRepMesh_SelfIntersectingWire;  return;
 406           }
 407           //look for intersection of two linear segments
 408           if(Abs(AB) <= RESOLUTION) continue;  //current segments seem parallel - no intersection
 409           //calculate coordinates of point of the intersection
 410           BC = B1*C2 - B2*C1;  xc = BC/AB;
 411           CA = C1*A2 - C2*A1;  yc = CA/AB;
 412           //check on belonging of intersection point to the both of segments
 413           if( Abs(xc-x11) > RESOLUTION && Abs(xc-x12) > RESOLUTION &&
 414               Abs(yc-y11) > RESOLUTION && Abs(yc-y12) > RESOLUTION &&
 415               Abs(xc-x21) > RESOLUTION && Abs(xc-x22) > RESOLUTION &&
 416               Abs(yc-y21) > RESOLUTION && Abs(yc-y22) > RESOLUTION )
 417           {
 418             if((xc-x11)*(xc-x12) < 0.0 && (yc-y11)*(yc-y12) < 0.0 &&
 419                (xc-x21)*(xc-x22) < 0.0 && (yc-y21)*(yc-y22) < 0.0)
 420             {
 421               //different segments may have common vertex (why "<" but "<=")
 422               myState = BRepMesh_SelfIntersectingWire;  return;
 423             }
 424           }
 425         }
 426       }
 427     }
 428   }
 429
 430   // Find holes
 431   for (i = nbwires; i >= 1; i--)
 432   {
 433     NbBiPoint = aWirePoints.Length() - aWireLength(i) + 1;
 434     aWirePoints.Split(NbBiPoint, aWire);
 435     AnalizeWire(aWire, Umin, Umax, Vmin, Vmax);
 436   }
 437 }
 438
 439
 440 //=======================================================================
 441 //function : Perform
 442 //purpose  :
 443 //=======================================================================
 444
 445 TopAbs_State BRepMesh_Classifier::Perform(const gp_Pnt2d& aPoint) const
 446 {
 447   Standard_Boolean isOut = Standard_False;
 448
 449   Standard_Integer cur, i, nb = TabClass.Length();
 450
 451   for (i = 1; i <= nb; i++)
 452   {
 453     cur = ((CSLib_Class2d*)TabClass(i))->SiDans(aPoint);
 454     if (cur == 0)
 455     {
 456       // Point is ON, but mark it as OUT
 457       isOut = Standard_True;
 458     }
 459     else
 460     {
 461       isOut = TabOrien(i)? (cur == -1) : (cur == 1);
 462     }
 463     if (isOut) return TopAbs_OUT;
 464   }
 465
 466   return TopAbs_IN;
 467 }
 468
 469
 470 //=======================================================================
 471 //function : Destroy
 472 //purpose  :
 473 //=======================================================================
 474
 475 void BRepMesh_Classifier::Destroy()
 476 {
 477   Standard_Integer i, nb = TabClass.Length();
 478   for (i = 1; i <= nb; i++)
 479   {
 480     if (TabClass(i))
 481     {
 482       delete ((CSLib_Class2d*)TabClass(i));
 483       TabClass(i) = NULL;
 484     }
 485   }
 486 }