2c64b046b16cca2421daa571c6d6d3a8fd28d9b5
[occt.git] / src / BRepMesh / BRepMesh_Classifier.cxx
1 // Created on: 1997-06-26
2 // Created by: Laurent PAINNOT
3 // Copyright (c) 1997-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #include <BRepMesh_Classifier.ixx>
18
19 // Kernel
20 #include <Precision.hxx>
21 #include <Standard_ErrorHandler.hxx>
22 #include <TColStd_ListOfTransient.hxx>
23 #include <TColStd_Array1OfInteger.hxx>
24 #include <TColStd_DataMapOfIntegerInteger.hxx>
25 #include <ElCLib.hxx>
26 // Geometry
27 #include <gp_Pnt.hxx>
28 #include <gp_Pnt2d.hxx>
29 #include <TColgp_SequenceOfPnt2d.hxx>
30 #include <TColgp_Array1OfPnt2d.hxx>
31 #include <GeomAbs_SurfaceType.hxx>
32 #include <Geom2dInt_Geom2dCurveTool.hxx>
33 #include <Geom2d_Line.hxx>
34 #include <Geom2d_BezierCurve.hxx>
35 #include <Geom2d_BSplineCurve.hxx>
36 #include <Geom2d_TrimmedCurve.hxx>
37 // Topology
38 #include <BRep_Tool.hxx>
39 #include <BRepTools.hxx>
40 #include <BRepTools_WireExplorer.hxx>
41 #include <BRepAdaptor_Curve2d.hxx>
42 #include <TopAbs_Orientation.hxx>
43 #include <TopExp.hxx>
44 #include <TopExp_Explorer.hxx>
45 #include <TopoDS_Edge.hxx>
46 #include <TopoDS.hxx>
47 #include <CSLib_Class2d.hxx>
48 #include <Poly_PolygonOnTriangulation.hxx>
49 // BRepMesh
50 #include <BRepMesh_Vertex.hxx>
51 #include <BRepMesh_Array1OfBiPoint.hxx>
52 #include <BRepMesh_PairOfPolygon.hxx>
53
54 static const Standard_Real PARALL_COND = Sin(M_PI/3.0);
55 static const Standard_Real RESOLUTION = 1.0E-16;
56
57 // Real mesh is created in the grid 10E5x10E5, so intersection
58 // should be cheched with double of discretization.
59 static const Standard_Real MIN_DIST = 2.E-5;
60
61 //=======================================================================
62 //function : AnalizeWire
63 //purpose  : 
64 //=======================================================================
65 void BRepMesh_Classifier::AnalizeWire (const TColgp_SequenceOfPnt2d&  theSeqPnt2d,
66                                        const Standard_Real theUmin,  const Standard_Real theUmax,
67                                        const Standard_Real theVmin,  const Standard_Real theVmax)
68 {
69   const Standard_Integer aNbPnts = theSeqPnt2d.Length();
70   if (aNbPnts < 2)
71     return;
72
73   // Accumulate angle
74   TColgp_Array1OfPnt2d aPClass(1, aNbPnts);
75   Standard_Real anAngle = 0.0;
76   gp_Pnt2d p1 = theSeqPnt2d(1), p2 = theSeqPnt2d(2), p3;
77   aPClass(1) = p1;
78   aPClass(2) = p2;
79   for (Standard_Integer i = 1; i <= aNbPnts; i++)
80   { 
81     Standard_Integer ii = i + 2;
82     if (ii > aNbPnts)
83     {
84       p3 = aPClass(ii - aNbPnts);
85     }
86     else
87     {
88       p3 = theSeqPnt2d.Value(ii);
89       aPClass(ii) = p3;
90     }
91
92     gp_Vec2d A(p1,p2), B(p2,p3);
93     if (A.SquareMagnitude() > 1.e-16 && B.SquareMagnitude() > 1.e-16)
94     {
95       const Standard_Real aCurAngle    = A.Angle(B);
96       const Standard_Real aCurAngleAbs = Abs(aCurAngle);
97       // Check if vectors are opposite
98       if (aCurAngleAbs > Precision::Angular() && (M_PI - aCurAngleAbs) > Precision::Angular())
99       {
100         anAngle += aCurAngle;
101         p1 = p2;
102       }
103     }
104     p2 = p3;
105   }
106   // Check for zero angle - treat self intersecting wire as outer
107   if (Abs(anAngle) < Precision::Angular())
108     anAngle = 0.0;
109
110   myTabClass.Append( (void *)new CSLib_Class2d(aPClass, myTolUV, myTolUV,
111                                                theUmin, theVmin, theUmax, theVmax) );
112   myTabOrient.Append( ((anAngle < 0.0) ? 0 : 1) );
113 }
114
115 //=======================================================================
116 //function : triangle2Area
117 //purpose  : calculating area under triangle
118 //=======================================================================
119 inline static Standard_Real triangle2Area(const gp_XY& p1, const gp_XY& p2)
120 {
121   return p1.Crossed(p2);
122 }
123
124 //=======================================================================
125 //function : getSegmentParams
126 //purpose  : extracting segment attributes 
127 //=======================================================================
128 static Standard_Real getSegmentParams(const BRepMesh_Array1OfBiPoint& theBiPoints,
129                                       const Standard_Integer Index,
130                                       Standard_Real& x11,
131                                       Standard_Real& y11,
132                                       Standard_Real& x12,
133                                       Standard_Real& y12,
134                                       Standard_Real& A,
135                                       Standard_Real& B,
136                                       Standard_Real& C)
137 {
138   Standard_Real *aCoordinates;
139   aCoordinates = ((Standard_Real*)(theBiPoints(Index).Coordinates()));
140   x11 =  aCoordinates[0];
141   y11 =  aCoordinates[1];
142   x12 =  aCoordinates[2];
143   y12 =  aCoordinates[3];
144   A   =  aCoordinates[5];
145   B   = -aCoordinates[4];
146   C   = - x11*A - y11*B;
147   return A*A+B*B;
148 }
149
150 //=======================================================================
151 //function : checkWiresIntersection
152 //purpose  : finding intersection.
153 //           If the intersection is found return Standard_True
154 //=======================================================================
155 static Standard_Boolean checkWiresIntersection(const Standard_Integer           theFirstWireId,
156                                                const Standard_Integer           theSecondWireId,
157                                                Standard_Integer* const          theFirstOuterSegmentId,
158                                                Standard_Integer                 theLastOuterSegmentId,
159                                                const TColStd_SequenceOfInteger& theWireLength,
160                                                const BRepMesh_Array1OfBiPoint&  theBiPoints,
161                                                const Standard_Boolean           findNextIntersection   = Standard_False,
162                                                const Standard_Boolean           isFirstSegment         = Standard_False,
163                                                Standard_Integer* const          theFirstInnerSegmentId = 0)
164 {
165   Standard_Real A1, B1, C1, A2, B2, C2, AB, BC, CA, xc, yc;
166   Standard_Real mu1, d, mu2;
167   Standard_Integer ik = *theFirstOuterSegmentId, jk;
168   Standard_Real x11, x12, y11, y12, x21, x22, y21, y22;
169
170   // Calculate bounds for first wire
171   Standard_Integer ikEnd = theLastOuterSegmentId;
172   Standard_Boolean isFirst = Standard_True;
173   if ( findNextIntersection )
174     isFirst = isFirstSegment;
175
176   // Calculate bounds for second wire
177   Standard_Integer jkStart = 0, jkEnd = 0;
178   for (jk = 1; jk <= theSecondWireId; jk++)
179   {
180     jkStart = jkEnd + 1;
181     jkEnd  += theWireLength(jk);
182   }
183
184   // total area under polygon (area of loop)
185   Standard_Real aLoopArea          = 0.0;
186   // area under first triangles of polygon
187   Standard_Real aFirstTriangleArea = 0.0;
188   // contains coordinates of the end point of segment if first intersection point is finding
189   // or coordinates of the intersecting point if second intersection point is finding
190   gp_XY aStartPoint;
191
192   for (; ik <= ikEnd; ik++)
193   {
194     mu1 = getSegmentParams(theBiPoints, ik, x11, y11, x12, y12, A1, B1, C1);
195     // for second intersection point we must count the area from first intersection point 
196     if ( !findNextIntersection )
197     {
198       aLoopArea = 0.0;
199       aStartPoint.SetCoord(x12, y12);
200     }
201
202     //for theFirstWireId == theSecondWireId the algorithm check current wire on selfintersection
203     if ( findNextIntersection && theFirstInnerSegmentId && isFirst)
204       jk = *theFirstInnerSegmentId;
205     else if (theSecondWireId == theFirstWireId)
206       jk = ik + 2;
207     else
208       jk = jkStart;
209
210     // Explore second wire
211     Standard_Boolean aFirstPass = Standard_True;
212     for (; jk <= jkEnd; jk++)
213     {
214       // don't check end's segment of the wire on selfrestriction
215       if ( theSecondWireId == theFirstWireId && isFirst && jk == ikEnd )
216         continue;
217
218       mu2 = getSegmentParams(theBiPoints, jk, x21, y21, x22, y22, A2, B2, C2);
219       gp_XY p2(x21, y21), p3(x22, y22);
220
221       //different segments may have common vertex (see OCC287 bug for example)
222       AB = A1*B2 - A2*B1;
223       //check on minimal of distance between current segment and points of another linear segments - OCC319
224       d = A1*x22 + B1*y22 + C1;
225       Standard_Real dTol = MIN_DIST*MIN_DIST;
226       if(theFirstWireId != theSecondWireId       && // if compared wires are different &&
227          AB*AB > PARALL_COND*PARALL_COND*mu1*mu2 && // angle between two segments greater then PARALL_COND &&
228          d*d   < dTol*mu1 &&                        // distance between vertex of the segment and other one's less then MIN_DIST
229          (x22-x11)*(x22-x12) < 0.0 && (y22-y11)*(y22-y12) < 0.0)
230       {
231         // if we finding the second intersection we must return Standard_False for setting
232         // self-intersection result flag
233         if ( findNextIntersection )
234           return Standard_False;
235
236         // we can step here when finding first intersection, return self-intersection flag
237         return Standard_True;
238       }
239
240       if( aFirstPass )
241         aFirstTriangleArea = triangle2Area(aStartPoint, p2);
242       
243       Standard_Real aTmpArea = triangle2Area(p2, p3);
244
245       //look for intersection of two linear segments
246       if(Abs(AB) <= RESOLUTION)
247       {
248         aLoopArea += aTmpArea;
249         continue;  //current segments seem parallel - no intersection
250       }
251       
252       //calculate coordinates of point of the intersection
253       BC = B1*C2 - B2*C1;  xc = BC/AB;
254       CA = C1*A2 - C2*A1;  yc = CA/AB;
255
256       // remember current intersection point and area of first triangle
257       if( findNextIntersection && ik == *theFirstOuterSegmentId && jk == *theFirstInnerSegmentId )
258       {
259         aStartPoint.SetCoord(xc, yc);
260         continue;
261       }
262
263       //check on belonging of intersection point to the both of segments
264       Standard_Boolean isOnLines = Standard_True;
265
266       Standard_Real dd[2][4] = { {(xc-x11), (xc-x12), (xc-x21), (xc-x22)},   //dX
267                                  {(yc-y11), (yc-y12), (yc-y21), (yc-y22)} }; //dY
268
269       for( Standard_Integer i = 0; i < 2; i++ )
270       {
271         if ( dd[i][0] * dd[i][1] > RESOLUTION || dd[i][2] * dd[i][3] > RESOLUTION )
272         {
273           isOnLines = Standard_False;
274           break;
275         }
276       }
277
278       // check the intersection point is on the ends of segments
279       if ( isOnLines )
280       {
281         for( Standard_Integer i = 0; i < 2; i++ )
282         {
283           //     if it's the last segment and intersection point lies at the end
284           if ( ( jk == jkEnd                                              ||
285           //     dX                        && dY
286           //     or when the start or the end point of the first segment
287                 (Abs(dd[0][0])  < MIN_DIST && Abs(dd[1][0])   < MIN_DIST) ||
288                 (Abs(dd[0][1])  < MIN_DIST && Abs(dd[1][1])   < MIN_DIST)) &&
289           //     is equal to one of the end points of the second
290                (Abs(dd[0][i+2]) < MIN_DIST && Abs(dd[1][i+2]) < MIN_DIST))
291           {
292             // no intersection
293             isOnLines = Standard_False;
294             aLoopArea = aTmpArea = 0.0;
295             aFirstPass = Standard_True;
296             break;
297           }
298         }
299       }
300
301
302       if( isOnLines )
303       {
304         p3.SetX(xc); p3.SetY(yc); 
305         aLoopArea += aFirstTriangleArea;             // First triangle area
306         aLoopArea += triangle2Area(p2, p3); 
307         aLoopArea += triangle2Area(p3, aStartPoint); // Last triangle area
308
309         if( Abs(aLoopArea)/2 > M_PI*MIN_DIST )
310         {
311           if ( findNextIntersection )
312           {
313             // intersection is found, but Standard_False returns, because area is too much
314             return Standard_False;
315           }
316
317           if ( checkWiresIntersection(theFirstWireId, theSecondWireId, &ik, ikEnd, theWireLength,
318                            theBiPoints, Standard_True, isFirst, &jk) )
319           {
320             // small crossing is not intersection, continue cheching
321             aLoopArea = aTmpArea = 0.0;
322             aFirstPass = Standard_True;
323           }
324           else
325           {
326             // if we found only one intersection
327             return Standard_True;
328           }
329         }
330         else if ( findNextIntersection )
331         {
332           // small intersection, skip double checking
333           *theFirstOuterSegmentId = ik;
334           *theFirstInnerSegmentId = jk + 1;
335           return Standard_True;
336         }
337       }
338       if ( aFirstPass )
339         aFirstPass = Standard_False;
340
341       aLoopArea += aTmpArea;
342     }
343     
344     if ( isFirst )
345       isFirst = Standard_False;
346   }
347   return Standard_False;
348 }
349
350
351 //=======================================================================
352 //function : BRepMesh_Classifier
353 //purpose  : 
354 //=======================================================================
355 BRepMesh_Classifier::BRepMesh_Classifier(const TopoDS_Face& theFace,
356                                          const Standard_Real theTolUV,
357                                          const BRepMesh_DataMapOfShapePairOfPolygon& theEdges,
358                                          const TColStd_IndexedMapOfInteger& theMap,
359                                          const Handle(BRepMesh_DataStructureOfDelaun)& theStructure,
360                                          const Standard_Real theUmin,
361                                          const Standard_Real theUmax,
362                                          const Standard_Real theVmin,
363                                          const Standard_Real theVmax)
364 : myTolUV( theTolUV ),
365   myFace ( theFace ),
366   myState( BRepMesh_NoError )
367 {
368   //-- impasse sur les surfs definies sur plus d une periode
369   //-- once definition
370   myFace.Orientation(TopAbs_FORWARD);
371   
372   TColgp_SequenceOfPnt2d    aWirePoints, aWire;
373   TColStd_SequenceOfInteger aWireLength;
374
375   TopoDS_Iterator aFaceExplorer;
376   for(aFaceExplorer.Initialize(myFace); aFaceExplorer.More(); aFaceExplorer.Next())
377   {
378     if(aFaceExplorer.Value().ShapeType() != TopAbs_WIRE)
379       continue;
380
381     // For each wire we create a data map, linking vertices (only
382     // the ends of edges) with their positions in the sequence of
383     // all 2d points from this wire.
384     // When we meet some vertex for the second time - the piece
385     // of sequence is treated for a HOLE and quits the sequence.
386     // Actually, we must unbind the vertices belonging to the
387     // loop from the map, but since they can't appear twice on the
388     // valid wire, leave them for a little speed up.
389     Standard_Integer aNbEdges = 0;
390     Standard_Integer aFirstIndex = 0, aLastIndex = 0;
391     Standard_Boolean isFalseWire = Standard_False;
392
393     TColgp_SequenceOfPnt2d aSeqPnt2d;
394     TColStd_DataMapOfIntegerInteger aNodeInSeq;
395
396     // Start traversing the wire
397     BRepTools_WireExplorer aWireExplorer;
398     for (aWireExplorer.Init(TopoDS::Wire( aFaceExplorer.Value() ), myFace); aWireExplorer.More(); aWireExplorer.Next())
399     {
400       TopoDS_Edge        anEdge   = aWireExplorer.Current();
401       TopAbs_Orientation anOrient = anEdge.Orientation();
402       if (anOrient != TopAbs_FORWARD && anOrient != TopAbs_REVERSED)
403         continue;
404
405       if (theEdges.IsBound(anEdge))
406       {
407         // Retrieve polygon
408         // Define the direction for adding points to aSeqPnt2d
409         Standard_Integer aIdxFirst, aIdxLast, aIdxIncr;
410
411         const BRepMesh_PairOfPolygon& aPair = theEdges.Find(anEdge);
412         Handle(Poly_PolygonOnTriangulation) aNOD;
413         if (anOrient == TopAbs_FORWARD)
414         {
415           aNOD = aPair.First();
416           aIdxFirst = 1;
417           aIdxLast  = aNOD->NbNodes();
418           aIdxIncr  = 1;
419         }
420         else
421         {
422           aNOD = aPair.Last();
423           aIdxFirst = aNOD->NbNodes();
424           aIdxLast  = 1;
425           aIdxIncr  = -1;
426         }
427         const TColStd_Array1OfInteger& anIndices = aNOD->Nodes();
428
429         // anIndexFirst and anIndexLast are the indices of first and last
430         // vertices of the edge in IndexedMap <Str>
431         const Standard_Integer anIndexFirst = theMap.FindKey( anIndices(aIdxFirst) );
432         const Standard_Integer anIndexLast  = theMap.FindKey( anIndices(aIdxLast) );
433
434         if (anIndexLast == anIndexFirst && (aIdxLast - aIdxFirst) == aIdxIncr)
435         {
436           // case of continuous set of degenerated edges
437           aLastIndex = anIndexLast;
438           continue;
439         }
440
441         // If there's a gap between edges -> raise <isFalseWire> flag
442         if (aNbEdges)
443         {
444           if (anIndexFirst != aLastIndex)
445           {
446             isFalseWire = Standard_True;
447             break;
448           }
449         }
450         else
451           aFirstIndex = anIndexFirst;
452
453         aLastIndex = anIndexLast;
454
455         // Record first vertex (to detect loops)
456         aNodeInSeq.Bind(anIndexFirst, (aSeqPnt2d.Length() + 1));
457
458         // Add vertices in sequence
459         for (Standard_Integer i = aIdxFirst; i != aIdxLast; i += aIdxIncr)
460         {
461           Standard_Integer anIndex = ((i == aIdxFirst) ? anIndexFirst : theMap.FindKey( anIndices(i) ));
462
463           gp_Pnt2d aPnt( theStructure->GetNode(anIndex).Coord() );
464           aSeqPnt2d.Append(aPnt);
465         }
466
467         // Now, is there a loop?
468         if (aNodeInSeq.IsBound(anIndexLast))
469         {
470           // Yes, treat it separately as a hole
471           // 1. Divide points into main wire and a loop
472           const Standard_Integer aIdxWireStart = aNodeInSeq(anIndexLast);
473           if(aIdxWireStart < aSeqPnt2d.Length())
474           {
475             aSeqPnt2d.Split(aIdxWireStart, aWire);
476             // 2. Proceed the loop
477             //AnalizeWire(aLoop, Umin, Umax, Vmin, Vmax, aWirePoints, aWireLength, NbBiPoint);
478             aWireLength.Append( aWire.Length() );
479             aWirePoints.Append( aWire );
480           }
481         }
482         aNbEdges++;
483       }
484     }
485
486     if (aNbEdges)
487     {
488       // Isn't it open?
489       if (isFalseWire || (aFirstIndex != aLastIndex) || aSeqPnt2d.Length() > 1)
490       {
491         myState = BRepMesh_OpenWire;
492         return;
493       }
494     }
495   }
496
497   const Standard_Integer aNbWires = aWireLength.Length();
498   Standard_Integer aNbBiPoint = aWirePoints.Length();
499   BRepMesh_Array1OfBiPoint aBiPoints(0, aNbBiPoint);
500   BRepMesh_BiPoint *aBiPoint = &(aBiPoints.ChangeValue(1));
501
502   // Fill array of segments (bi-points)
503   Standard_Integer k = 1;
504   for (Standard_Integer i = 1; i <= aNbWires; i++)
505   {
506     Standard_Real x1 = 0., y1 = 0., x2, y2, aXstart = 0., aYstart = 0.;
507     const Standard_Integer aLen = aWireLength(i) + 1;
508     for (Standard_Integer j = 1; j <= aLen; j++)
509     {
510       // Obtain last point of the segment
511       if (j == aLen)
512       {
513         x2 = aXstart;
514         y2 = aYstart;
515       }
516       else
517       {
518         const gp_Pnt2d& aPnt = aWirePoints(k);
519         k++;
520
521         x2 = aPnt.X();
522         y2 = aPnt.Y();
523       }
524       // Build segment (bi-point)
525       if (j == 1)
526       {
527         aXstart = x2;
528         aYstart = y2;
529       }
530       else
531       {
532         Standard_Real *aCoordinates1 = ((Standard_Real*)(aBiPoint->Coordinates()));
533         aBiPoint++;
534
535         aCoordinates1[0] = x1;
536         aCoordinates1[1] = y1;
537         aCoordinates1[2] = x2;
538         aCoordinates1[3] = y2;
539         aCoordinates1[4] = x2 - x1;
540         aCoordinates1[5] = y2 - y1;
541       }
542       x1 = x2;
543       y1 = y2;
544     }
545   }
546
547   // Search the intersection
548   // Explore first wire
549   Standard_Integer ikEnd = 0;
550   for(Standard_Integer i = 1; i <= aNbWires; i++)
551   {
552     Standard_Integer ik = ikEnd + 1;
553     ikEnd += aWireLength(i);
554
555     // Explore second wire
556     for (Standard_Integer j = i; j <= aNbWires; j++)
557     {
558       if ( checkWiresIntersection(i, j, &ik, ikEnd, aWireLength, aBiPoints) )
559       {
560         myState = BRepMesh_SelfIntersectingWire;
561         return;
562       }
563     }
564   }
565
566   // Find holes
567   for (Standard_Integer i = aNbWires; i >= 1; i--)
568   {
569     aNbBiPoint = aWirePoints.Length() - aWireLength(i) + 1;
570     aWirePoints.Split(aNbBiPoint, aWire);
571     AnalizeWire(aWire, theUmin, theUmax, theVmin, theVmax);
572   }
573 }
574
575
576 //=======================================================================
577 //function : Perform
578 //purpose  : 
579 //=======================================================================
580 TopAbs_State BRepMesh_Classifier::Perform(const gp_Pnt2d& thePoint) const
581 {
582   Standard_Boolean isOut = Standard_False;
583   Standard_Integer aNb   = myTabClass.Length();
584   
585   for (Standard_Integer i = 1; i <= aNb; i++)
586   {
587     Standard_Integer aCur = ((CSLib_Class2d*)myTabClass(i))->SiDans(thePoint);
588     if (aCur == 0)
589     {
590       // Point is ON, but mark it as OUT
591       isOut = Standard_True;
592     }
593     else
594       isOut = myTabOrient(i)? (aCur == -1) : (aCur == 1);
595     
596     if (isOut)
597       return TopAbs_OUT;
598   }
599
600   return TopAbs_IN;
601 }
602
603
604 //=======================================================================
605 //function : Destroy
606 //purpose  : 
607 //=======================================================================
608 void BRepMesh_Classifier::Destroy()
609 {
610   Standard_Integer aNb = myTabClass.Length();
611   for (Standard_Integer i = 1; i <= aNb; i++)
612   {
613     if (myTabClass(i))
614     {
615       delete ((CSLib_Class2d*)myTabClass(i));
616       myTabClass(i) = NULL;
617     }
618   }
619 }