1 // Created on: 1997-06-26
2 // Created by: Laurent PAINNOT
3 // Copyright (c) 1997-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
23 #include <BRepMesh_Classifier.ixx>
26 #include <Precision.hxx>
27 #include <Standard_ErrorHandler.hxx>
28 #include <TColStd_ListOfTransient.hxx>
29 #include <TColStd_Array1OfInteger.hxx>
30 #include <TColStd_DataMapOfIntegerInteger.hxx>
34 #include <gp_Pnt2d.hxx>
35 #include <TColgp_SequenceOfPnt2d.hxx>
36 #include <TColgp_Array1OfPnt2d.hxx>
37 #include <GeomAbs_SurfaceType.hxx>
38 #include <Geom2dInt_Geom2dCurveTool.hxx>
39 #include <Geom2d_Line.hxx>
40 #include <Geom2d_BezierCurve.hxx>
41 #include <Geom2d_BSplineCurve.hxx>
42 #include <Geom2d_TrimmedCurve.hxx>
44 #include <BRep_Tool.hxx>
45 #include <BRepTools.hxx>
46 #include <BRepTools_WireExplorer.hxx>
47 #include <BRepAdaptor_Curve2d.hxx>
48 #include <TopAbs_Orientation.hxx>
50 #include <TopExp_Explorer.hxx>
51 #include <TopoDS_Edge.hxx>
53 #include <CSLib_Class2d.hxx>
54 #include <Poly_PolygonOnTriangulation.hxx>
56 #include <BRepMesh_Vertex.hxx>
57 #include <BRepMesh_Array1OfBiPoint.hxx>
58 #include <BRepMesh_PairOfPolygon.hxx>
60 static const Standard_Real PARALL_COND = Sin(M_PI/3.0);
61 static const Standard_Real RESOLUTION = 1.0E-16;
63 // Real mesh is created in the grid 10E5x10E5, so intersection
64 // should be cheched with double of discretization.
65 static const Standard_Real MIN_DIST = 2.E-5;
67 //=======================================================================
70 //=======================================================================
71 static Standard_Boolean IsLine(const Handle(Geom2d_Curve)& theCurve2d)
73 Standard_Boolean IsALine = Standard_False;
74 if ( theCurve2d->IsKind( STANDARD_TYPE(Geom2d_Line) ) )
76 IsALine = Standard_True;
78 else if ( theCurve2d->IsKind( STANDARD_TYPE(Geom2d_BSplineCurve) ) )
80 Handle(Geom2d_BSplineCurve) aBSpline = *((Handle(Geom2d_BSplineCurve)*)&theCurve2d);
81 IsALine = (aBSpline->NbPoles() == 2);
83 else if ( theCurve2d->IsKind( STANDARD_TYPE(Geom2d_BezierCurve) ) )
85 Handle(Geom2d_BezierCurve) aBezier = *((Handle(Geom2d_BezierCurve)*)&theCurve2d);
86 IsALine = (aBezier->NbPoles() == 2);
88 else if ( theCurve2d->IsKind( STANDARD_TYPE(Geom2d_TrimmedCurve) ) )
90 Handle(Geom2d_TrimmedCurve) aTrimmedCurve = *((Handle(Geom2d_TrimmedCurve)*)&theCurve2d);
91 IsALine = IsLine(aTrimmedCurve->BasisCurve());
96 //=======================================================================
97 //function : AnalizeWire
99 //=======================================================================
100 void BRepMesh_Classifier::AnalizeWire (const TColgp_SequenceOfPnt2d& theSeqPnt2d,
101 const Standard_Real theUmin, const Standard_Real theUmax,
102 const Standard_Real theVmin, const Standard_Real theVmax)
104 const Standard_Integer aNbPnts = theSeqPnt2d.Length();
109 TColgp_Array1OfPnt2d aPClass(1, aNbPnts);
110 Standard_Real anAngle = 0.0;
111 gp_Pnt2d p1 = theSeqPnt2d(1), p2 = theSeqPnt2d(2), p3;
114 for (Standard_Integer i = 1; i <= aNbPnts; i++)
116 Standard_Integer ii = i + 2;
119 p3 = aPClass(ii - aNbPnts);
123 p3 = theSeqPnt2d.Value(ii);
127 gp_Vec2d A(p1,p2), B(p2,p3);
128 if (A.SquareMagnitude() > 1.e-16 && B.SquareMagnitude() > 1.e-16)
130 const Standard_Real aCurAngle = A.Angle(B);
131 const Standard_Real aCurAngleAbs = Abs(aCurAngle);
132 // Check if vectors are opposite
133 if (aCurAngleAbs > Precision::Angular() && (M_PI - aCurAngleAbs) > Precision::Angular())
135 anAngle += aCurAngle;
141 // Check for zero angle - treat self intersecting wire as outer
142 if (Abs(anAngle) < Precision::Angular())
145 myTabClass.Append( (void *)new CSLib_Class2d(aPClass, myTolUV, myTolUV,
146 theUmin, theVmin, theUmax, theVmax) );
147 myTabOrient.Append( ((anAngle < 0.0) ? 0 : 1) );
150 //=======================================================================
151 //function : triangle2Area
152 //purpose : calculating area under triangle
153 //=======================================================================
154 inline static Standard_Real triangle2Area(const gp_XY& p1, const gp_XY& p2)
156 return p1.Crossed(p2);
159 //=======================================================================
160 //function : getSegmentParams
161 //purpose : extracting segment attributes
162 //=======================================================================
163 static Standard_Real getSegmentParams(const BRepMesh_Array1OfBiPoint& theBiPoints,
164 const Standard_Integer Index,
173 Standard_Real *aCoordinates;
174 aCoordinates = ((Standard_Real*)(theBiPoints(Index).Coordinates()));
175 x11 = aCoordinates[0];
176 y11 = aCoordinates[1];
177 x12 = aCoordinates[2];
178 y12 = aCoordinates[3];
180 B = -aCoordinates[4];
185 //=======================================================================
186 //function : checkWiresIntersection
187 //purpose : finding intersection.
188 // If the intersection is found return Standard_True
189 //=======================================================================
190 static Standard_Boolean checkWiresIntersection(const Standard_Integer theFirstWireId,
191 const Standard_Integer theSecondWireId,
192 Standard_Integer* const theFirstOuterSegmentId,
193 Standard_Integer theLastOuterSegmentId,
194 const TColStd_SequenceOfInteger& theWireLength,
195 const BRepMesh_Array1OfBiPoint& theBiPoints,
196 const Standard_Boolean findNextIntersection = Standard_False,
197 const Standard_Boolean isFirstSegment = Standard_False,
198 Standard_Integer* const theFirstInnerSegmentId = 0)
200 Standard_Real A1, B1, C1, A2, B2, C2, AB, BC, CA, xc, yc;
201 Standard_Real mu1, d, mu2;
202 Standard_Integer ik = *theFirstOuterSegmentId, jk;
203 Standard_Real x11, x12, y11, y12, x21, x22, y21, y22;
205 // Calculate bounds for first wire
206 Standard_Integer ikEnd = theLastOuterSegmentId;
207 Standard_Boolean isFirst = Standard_True;
208 if ( findNextIntersection )
209 isFirst = isFirstSegment;
211 // Calculate bounds for second wire
212 Standard_Integer jkStart = 0, jkEnd = 0;
213 for (jk = 1; jk <= theSecondWireId; jk++)
216 jkEnd += theWireLength(jk);
219 // total area under polygon (area of loop)
220 Standard_Real aLoopArea = 0.0;
221 // area under first triangles of polygon
222 Standard_Real aFirstTriangleArea = 0.0;
223 // contains coordinates of the end point of segment if first intersection point is finding
224 // or coordinates of the intersecting point if second intersection point is finding
227 for (; ik <= ikEnd; ik++)
229 mu1 = getSegmentParams(theBiPoints, ik, x11, y11, x12, y12, A1, B1, C1);
230 // for second intersection point we must count the area from first intersection point
231 if ( !findNextIntersection )
234 aStartPoint.SetCoord(x12, y12);
237 //for theFirstWireId == theSecondWireId the algorithm check current wire on selfintersection
238 if ( findNextIntersection && theFirstInnerSegmentId && isFirst)
239 jk = *theFirstInnerSegmentId;
240 else if (theSecondWireId == theFirstWireId)
245 // Explore second wire
246 Standard_Boolean aFirstPass = Standard_True;
247 for (; jk <= jkEnd; jk++)
249 // don't check end's segment of the wire on selfrestriction
250 if ( theSecondWireId == theFirstWireId && isFirst && jk == ikEnd )
253 mu2 = getSegmentParams(theBiPoints, jk, x21, y21, x22, y22, A2, B2, C2);
254 gp_XY p2(x21, y21), p3(x22, y22);
256 //different segments may have common vertex (see OCC287 bug for example)
258 //check on minimal of distance between current segment and points of another linear segments - OCC319
259 d = A1*x22 + B1*y22 + C1;
260 Standard_Real dTol = MIN_DIST*MIN_DIST;
261 if(theFirstWireId != theSecondWireId && // if compared wires are different &&
262 AB*AB > PARALL_COND*PARALL_COND*mu1*mu2 && // angle between two segments greater then PARALL_COND &&
263 d*d < dTol*mu1 && // distance between vertex of the segment and other one's less then MIN_DIST
264 (x22-x11)*(x22-x12) < 0.0 && (y22-y11)*(y22-y12) < 0.0)
266 // if we finding the second intersection we must return Standard_False for setting
267 // self-intersection result flag
268 if ( findNextIntersection )
269 return Standard_False;
271 // we can step here when finding first intersection, return self-intersection flag
272 return Standard_True;
276 aFirstTriangleArea = triangle2Area(aStartPoint, p2);
278 Standard_Real aTmpArea = triangle2Area(p2, p3);
280 //look for intersection of two linear segments
281 if(Abs(AB) <= RESOLUTION)
283 aLoopArea += aTmpArea;
284 continue; //current segments seem parallel - no intersection
287 //calculate coordinates of point of the intersection
288 BC = B1*C2 - B2*C1; xc = BC/AB;
289 CA = C1*A2 - C2*A1; yc = CA/AB;
291 // remember current intersection point and area of first triangle
292 if( findNextIntersection && ik == *theFirstOuterSegmentId && jk == *theFirstInnerSegmentId )
294 aStartPoint.SetCoord(xc, yc);
298 //check on belonging of intersection point to the both of segments
299 Standard_Boolean isOnLines = Standard_True;
301 Standard_Real dd[2][4] = { {(xc-x11), (xc-x12), (xc-x21), (xc-x22)}, //dX
302 {(yc-y11), (yc-y12), (yc-y21), (yc-y22)} }; //dY
304 Standard_Integer i = 0;
307 if ( dd[i][0]*dd[i][1] > dTol || dd[i][2]*dd[i][3] > dTol)
309 isOnLines = Standard_False;
314 // check the intersection point is on the ends of segments
317 for( i = 0; i < 2; i++ )
319 // if it's the last segment and intersection point lies at the end
320 if ( ( jk == jkEnd ||
322 // or when the start or the end point of the first segment
323 (Abs(dd[0][0]) < MIN_DIST && Abs(dd[1][0]) < MIN_DIST) ||
324 (Abs(dd[0][1]) < MIN_DIST && Abs(dd[1][1]) < MIN_DIST)) &&
325 // is equal to one of the end points of the second
326 (Abs(dd[0][i+2]) < MIN_DIST && Abs(dd[1][i+2]) < MIN_DIST))
329 isOnLines = Standard_False;
330 aLoopArea = aTmpArea = 0.0;
331 aFirstPass = Standard_True;
340 p3.SetX(xc); p3.SetY(yc);
341 aLoopArea += aFirstTriangleArea; // First triangle area
342 aLoopArea += triangle2Area(p2, p3);
343 aLoopArea += triangle2Area(p3, aStartPoint); // Last triangle area
345 if( Abs(aLoopArea)/2 > M_PI*MIN_DIST )
347 if ( findNextIntersection )
349 // intersection is found, but Standard_False returns, because area is too much
350 return Standard_False;
353 if ( checkWiresIntersection(theFirstWireId, theSecondWireId, &ik, ikEnd, theWireLength,
354 theBiPoints, Standard_True, isFirst, &jk) )
356 // small crossing is not intersection, continue cheching
357 aLoopArea = aTmpArea = 0.0;
358 aFirstPass = Standard_True;
362 // if we found only one intersection
363 return Standard_True;
366 else if ( findNextIntersection )
368 // small intersection, skip double checking
369 *theFirstOuterSegmentId = ik;
370 *theFirstInnerSegmentId = jk + 1;
371 return Standard_True;
375 aFirstPass = Standard_False;
377 aLoopArea += aTmpArea;
381 isFirst = Standard_False;
383 return Standard_False;
387 //=======================================================================
388 //function : BRepMesh_Classifier
390 //=======================================================================
391 BRepMesh_Classifier::BRepMesh_Classifier(const TopoDS_Face& theFace,
392 const Standard_Real theTolUV,
393 const BRepMesh_DataMapOfShapePairOfPolygon& theEdges,
394 const TColStd_IndexedMapOfInteger& theMap,
395 const Handle(BRepMesh_DataStructureOfDelaun)& theStructure,
396 const Standard_Real theUmin,
397 const Standard_Real theUmax,
398 const Standard_Real theVmin,
399 const Standard_Real theVmax)
400 : myTolUV( theTolUV ),
402 myState( BRepMesh_NoError )
404 //-- impasse sur les surfs definies sur plus d une periode
406 myFace.Orientation(TopAbs_FORWARD);
408 TColgp_SequenceOfPnt2d aWirePoints, aWire;
409 TColStd_SequenceOfInteger aWireLength;
411 TopoDS_Iterator aFaceExplorer;
412 for(aFaceExplorer.Initialize(myFace); aFaceExplorer.More(); aFaceExplorer.Next())
414 if(aFaceExplorer.Value().ShapeType() != TopAbs_WIRE)
417 // For each wire we create a data map, linking vertices (only
418 // the ends of edges) with their positions in the sequence of
419 // all 2d points from this wire.
420 // When we meet some vertex for the second time - the piece
421 // of sequence is treated for a HOLE and quits the sequence.
422 // Actually, we must unbind the vertices belonging to the
423 // loop from the map, but since they can't appear twice on the
424 // valid wire, leave them for a little speed up.
425 Standard_Integer aNbEdges = 0;
426 Standard_Integer aFirstIndex = 0, aLastIndex = 0;
427 Standard_Boolean isFalseWire = Standard_False;
429 TColgp_SequenceOfPnt2d aSeqPnt2d;
430 TColStd_DataMapOfIntegerInteger aNodeInSeq;
432 // Start traversing the wire
433 BRepTools_WireExplorer aWireExplorer;
434 for (aWireExplorer.Init(TopoDS::Wire( aFaceExplorer.Value() ), myFace); aWireExplorer.More(); aWireExplorer.Next())
436 TopoDS_Edge anEdge = aWireExplorer.Current();
437 TopAbs_Orientation anOrient = anEdge.Orientation();
438 if (anOrient != TopAbs_FORWARD && anOrient != TopAbs_REVERSED)
441 if (theEdges.IsBound(anEdge))
444 // Define the direction for adding points to aSeqPnt2d
445 Standard_Integer aIdxFirst, aIdxLast, aIdxIncr;
447 const BRepMesh_PairOfPolygon& aPair = theEdges.Find(anEdge);
448 Handle(Poly_PolygonOnTriangulation) aNOD;
449 if (anOrient == TopAbs_FORWARD)
451 aNOD = aPair.First();
453 aIdxLast = aNOD->NbNodes();
459 aIdxFirst = aNOD->NbNodes();
463 const TColStd_Array1OfInteger& anIndices = aNOD->Nodes();
465 // anIndexFirst and anIndexLast are the indices of first and last
466 // vertices of the edge in IndexedMap <Str>
467 const Standard_Integer anIndexFirst = theMap.FindKey( anIndices(aIdxFirst) );
468 const Standard_Integer anIndexLast = theMap.FindKey( anIndices(aIdxLast) );
470 if (anIndexLast == anIndexFirst && (aIdxLast - aIdxFirst) == aIdxIncr)
472 // case of continuous set of degenerated edges
473 aLastIndex = anIndexLast;
477 // If there's a gap between edges -> raise <isFalseWire> flag
480 if (anIndexFirst != aLastIndex)
482 isFalseWire = Standard_True;
487 aFirstIndex = anIndexFirst;
489 aLastIndex = anIndexLast;
491 // Record first vertex (to detect loops)
492 aNodeInSeq.Bind(anIndexFirst, (aSeqPnt2d.Length() + 1));
494 // Add vertices in sequence
495 for (Standard_Integer i = aIdxFirst; i != aIdxLast; i += aIdxIncr)
497 Standard_Integer anIndex = ((i == aIdxFirst) ? anIndexFirst : theMap.FindKey( anIndices(i) ));
499 gp_Pnt2d aPnt( theStructure->GetNode(anIndex).Coord() );
500 aSeqPnt2d.Append(aPnt);
503 // Now, is there a loop?
504 if (aNodeInSeq.IsBound(anIndexLast))
506 // Yes, treat it separately as a hole
507 // 1. Divide points into main wire and a loop
508 const Standard_Integer aIdxWireStart = aNodeInSeq(anIndexLast);
509 if(aIdxWireStart < aSeqPnt2d.Length())
511 aSeqPnt2d.Split(aIdxWireStart, aWire);
512 // 2. Proceed the loop
513 //AnalizeWire(aLoop, Umin, Umax, Vmin, Vmax, aWirePoints, aWireLength, NbBiPoint);
514 aWireLength.Append( aWire.Length() );
515 aWirePoints.Append( aWire );
525 if (isFalseWire || (aFirstIndex != aLastIndex) || aSeqPnt2d.Length() > 1)
527 myState = BRepMesh_OpenWire;
533 const Standard_Integer aNbWires = aWireLength.Length();
534 Standard_Integer aNbBiPoint = aWirePoints.Length();
535 BRepMesh_Array1OfBiPoint aBiPoints(0, aNbBiPoint);
536 BRepMesh_BiPoint *aBiPoint = &(aBiPoints.ChangeValue(1));
538 // Fill array of segments (bi-points)
539 Standard_Integer k = 1;
540 for (Standard_Integer i = 1; i <= aNbWires; i++)
542 Standard_Real x1, y1, x2, y2, aXstart, aYstart;
543 const Standard_Integer aLen = aWireLength(i) + 1;
544 for (Standard_Integer j = 1; j <= aLen; j++)
546 // Obtain last point of the segment
554 const gp_Pnt2d& aPnt = aWirePoints(k);
560 // Build segment (bi-point)
568 Standard_Real *aCoordinates1 = ((Standard_Real*)(aBiPoint->Coordinates()));
571 aCoordinates1[0] = x1;
572 aCoordinates1[1] = y1;
573 aCoordinates1[2] = x2;
574 aCoordinates1[3] = y2;
575 aCoordinates1[4] = x2 - x1;
576 aCoordinates1[5] = y2 - y1;
583 // Search the intersection
584 // Explore first wire
585 Standard_Integer ikEnd = 0;
586 for(Standard_Integer i = 1; i <= aNbWires; i++)
588 Standard_Integer ik = ikEnd + 1;
589 ikEnd += aWireLength(i);
591 // Explore second wire
592 for (Standard_Integer j = i; j <= aNbWires; j++)
594 if ( checkWiresIntersection(i, j, &ik, ikEnd, aWireLength, aBiPoints) )
596 myState = BRepMesh_SelfIntersectingWire;
603 for (Standard_Integer i = aNbWires; i >= 1; i--)
605 aNbBiPoint = aWirePoints.Length() - aWireLength(i) + 1;
606 aWirePoints.Split(aNbBiPoint, aWire);
607 AnalizeWire(aWire, theUmin, theUmax, theVmin, theVmax);
612 //=======================================================================
615 //=======================================================================
616 TopAbs_State BRepMesh_Classifier::Perform(const gp_Pnt2d& thePoint) const
618 Standard_Boolean isOut = Standard_False;
619 Standard_Integer aNb = myTabClass.Length();
621 for (Standard_Integer i = 1; i <= aNb; i++)
623 Standard_Integer aCur = ((CSLib_Class2d*)myTabClass(i))->SiDans(thePoint);
626 // Point is ON, but mark it as OUT
627 isOut = Standard_True;
630 isOut = myTabOrient(i)? (aCur == -1) : (aCur == 1);
640 //=======================================================================
643 //=======================================================================
644 void BRepMesh_Classifier::Destroy()
646 Standard_Integer aNb = myTabClass.Length();
647 for (Standard_Integer i = 1; i <= aNb; i++)
651 delete ((CSLib_Class2d*)myTabClass(i));
652 myTabClass(i) = NULL;