0024428: Implementation of LGPL license
[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
9 // under the terms of the GNU Lesser General Public 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       Standard_Integer i = 0;
270       for(; i < 2; i++ )
271       {
272         if ( dd[i][0]*dd[i][1] > dTol || dd[i][2]*dd[i][3] > dTol)
273         {
274           isOnLines = Standard_False;
275           break;
276         }
277       }
278
279       // check the intersection point is on the ends of segments
280       if ( isOnLines )
281       {
282         for( i = 0; i < 2; i++ )
283         {
284           //     if it's the last segment and intersection point lies at the end
285           if ( ( jk == jkEnd                                              ||
286           //     dX                        && dY
287           //     or when the start or the end point of the first segment
288                 (Abs(dd[0][0])  < MIN_DIST && Abs(dd[1][0])   < MIN_DIST) ||
289                 (Abs(dd[0][1])  < MIN_DIST && Abs(dd[1][1])   < MIN_DIST)) &&
290           //     is equal to one of the end points of the second
291                (Abs(dd[0][i+2]) < MIN_DIST && Abs(dd[1][i+2]) < MIN_DIST))
292           {
293             // no intersection
294             isOnLines = Standard_False;
295             aLoopArea = aTmpArea = 0.0;
296             aFirstPass = Standard_True;
297             break;
298           }
299         }
300       }
301
302
303       if( isOnLines )
304       {
305         p3.SetX(xc); p3.SetY(yc); 
306         aLoopArea += aFirstTriangleArea;             // First triangle area
307         aLoopArea += triangle2Area(p2, p3); 
308         aLoopArea += triangle2Area(p3, aStartPoint); // Last triangle area
309
310         if( Abs(aLoopArea)/2 > M_PI*MIN_DIST )
311         {
312           if ( findNextIntersection )
313           {
314             // intersection is found, but Standard_False returns, because area is too much
315             return Standard_False;
316           }
317
318           if ( checkWiresIntersection(theFirstWireId, theSecondWireId, &ik, ikEnd, theWireLength,
319                            theBiPoints, Standard_True, isFirst, &jk) )
320           {
321             // small crossing is not intersection, continue cheching
322             aLoopArea = aTmpArea = 0.0;
323             aFirstPass = Standard_True;
324           }
325           else
326           {
327             // if we found only one intersection
328             return Standard_True;
329           }
330         }
331         else if ( findNextIntersection )
332         {
333           // small intersection, skip double checking
334           *theFirstOuterSegmentId = ik;
335           *theFirstInnerSegmentId = jk + 1;
336           return Standard_True;
337         }
338       }
339       if ( aFirstPass )
340         aFirstPass = Standard_False;
341
342       aLoopArea += aTmpArea;
343     }
344     
345     if ( isFirst )
346       isFirst = Standard_False;
347   }
348   return Standard_False;
349 }
350
351
352 //=======================================================================
353 //function : BRepMesh_Classifier
354 //purpose  : 
355 //=======================================================================
356 BRepMesh_Classifier::BRepMesh_Classifier(const TopoDS_Face& theFace,
357                                          const Standard_Real theTolUV,
358                                          const BRepMesh_DataMapOfShapePairOfPolygon& theEdges,
359                                          const TColStd_IndexedMapOfInteger& theMap,
360                                          const Handle(BRepMesh_DataStructureOfDelaun)& theStructure,
361                                          const Standard_Real theUmin,
362                                          const Standard_Real theUmax,
363                                          const Standard_Real theVmin,
364                                          const Standard_Real theVmax)
365 : myTolUV( theTolUV ),
366   myFace ( theFace ),
367   myState( BRepMesh_NoError )
368 {
369   //-- impasse sur les surfs definies sur plus d une periode
370   //-- once definition
371   myFace.Orientation(TopAbs_FORWARD);
372   
373   TColgp_SequenceOfPnt2d    aWirePoints, aWire;
374   TColStd_SequenceOfInteger aWireLength;
375
376   TopoDS_Iterator aFaceExplorer;
377   for(aFaceExplorer.Initialize(myFace); aFaceExplorer.More(); aFaceExplorer.Next())
378   {
379     if(aFaceExplorer.Value().ShapeType() != TopAbs_WIRE)
380       continue;
381
382     // For each wire we create a data map, linking vertices (only
383     // the ends of edges) with their positions in the sequence of
384     // all 2d points from this wire.
385     // When we meet some vertex for the second time - the piece
386     // of sequence is treated for a HOLE and quits the sequence.
387     // Actually, we must unbind the vertices belonging to the
388     // loop from the map, but since they can't appear twice on the
389     // valid wire, leave them for a little speed up.
390     Standard_Integer aNbEdges = 0;
391     Standard_Integer aFirstIndex = 0, aLastIndex = 0;
392     Standard_Boolean isFalseWire = Standard_False;
393
394     TColgp_SequenceOfPnt2d aSeqPnt2d;
395     TColStd_DataMapOfIntegerInteger aNodeInSeq;
396
397     // Start traversing the wire
398     BRepTools_WireExplorer aWireExplorer;
399     for (aWireExplorer.Init(TopoDS::Wire( aFaceExplorer.Value() ), myFace); aWireExplorer.More(); aWireExplorer.Next())
400     {
401       TopoDS_Edge        anEdge   = aWireExplorer.Current();
402       TopAbs_Orientation anOrient = anEdge.Orientation();
403       if (anOrient != TopAbs_FORWARD && anOrient != TopAbs_REVERSED)
404         continue;
405
406       if (theEdges.IsBound(anEdge))
407       {
408         // Retrieve polygon
409         // Define the direction for adding points to aSeqPnt2d
410         Standard_Integer aIdxFirst, aIdxLast, aIdxIncr;
411
412         const BRepMesh_PairOfPolygon& aPair = theEdges.Find(anEdge);
413         Handle(Poly_PolygonOnTriangulation) aNOD;
414         if (anOrient == TopAbs_FORWARD)
415         {
416           aNOD = aPair.First();
417           aIdxFirst = 1;
418           aIdxLast  = aNOD->NbNodes();
419           aIdxIncr  = 1;
420         }
421         else
422         {
423           aNOD = aPair.Last();
424           aIdxFirst = aNOD->NbNodes();
425           aIdxLast  = 1;
426           aIdxIncr  = -1;
427         }
428         const TColStd_Array1OfInteger& anIndices = aNOD->Nodes();
429
430         // anIndexFirst and anIndexLast are the indices of first and last
431         // vertices of the edge in IndexedMap <Str>
432         const Standard_Integer anIndexFirst = theMap.FindKey( anIndices(aIdxFirst) );
433         const Standard_Integer anIndexLast  = theMap.FindKey( anIndices(aIdxLast) );
434
435         if (anIndexLast == anIndexFirst && (aIdxLast - aIdxFirst) == aIdxIncr)
436         {
437           // case of continuous set of degenerated edges
438           aLastIndex = anIndexLast;
439           continue;
440         }
441
442         // If there's a gap between edges -> raise <isFalseWire> flag
443         if (aNbEdges)
444         {
445           if (anIndexFirst != aLastIndex)
446           {
447             isFalseWire = Standard_True;
448             break;
449           }
450         }
451         else
452           aFirstIndex = anIndexFirst;
453
454         aLastIndex = anIndexLast;
455
456         // Record first vertex (to detect loops)
457         aNodeInSeq.Bind(anIndexFirst, (aSeqPnt2d.Length() + 1));
458
459         // Add vertices in sequence
460         for (Standard_Integer i = aIdxFirst; i != aIdxLast; i += aIdxIncr)
461         {
462           Standard_Integer anIndex = ((i == aIdxFirst) ? anIndexFirst : theMap.FindKey( anIndices(i) ));
463
464           gp_Pnt2d aPnt( theStructure->GetNode(anIndex).Coord() );
465           aSeqPnt2d.Append(aPnt);
466         }
467
468         // Now, is there a loop?
469         if (aNodeInSeq.IsBound(anIndexLast))
470         {
471           // Yes, treat it separately as a hole
472           // 1. Divide points into main wire and a loop
473           const Standard_Integer aIdxWireStart = aNodeInSeq(anIndexLast);
474           if(aIdxWireStart < aSeqPnt2d.Length())
475           {
476             aSeqPnt2d.Split(aIdxWireStart, aWire);
477             // 2. Proceed the loop
478             //AnalizeWire(aLoop, Umin, Umax, Vmin, Vmax, aWirePoints, aWireLength, NbBiPoint);
479             aWireLength.Append( aWire.Length() );
480             aWirePoints.Append( aWire );
481           }
482         }
483         aNbEdges++;
484       }
485     }
486
487     if (aNbEdges)
488     {
489       // Isn't it open?
490       if (isFalseWire || (aFirstIndex != aLastIndex) || aSeqPnt2d.Length() > 1)
491       {
492         myState = BRepMesh_OpenWire;
493         return;
494       }
495     }
496   }
497
498   const Standard_Integer aNbWires = aWireLength.Length();
499   Standard_Integer aNbBiPoint = aWirePoints.Length();
500   BRepMesh_Array1OfBiPoint aBiPoints(0, aNbBiPoint);
501   BRepMesh_BiPoint *aBiPoint = &(aBiPoints.ChangeValue(1));
502
503   // Fill array of segments (bi-points)
504   Standard_Integer k = 1;
505   for (Standard_Integer i = 1; i <= aNbWires; i++)
506   {
507     Standard_Real x1 = 0., y1 = 0., x2, y2, aXstart = 0., aYstart = 0.;
508     const Standard_Integer aLen = aWireLength(i) + 1;
509     for (Standard_Integer j = 1; j <= aLen; j++)
510     {
511       // Obtain last point of the segment
512       if (j == aLen)
513       {
514         x2 = aXstart;
515         y2 = aYstart;
516       }
517       else
518       {
519         const gp_Pnt2d& aPnt = aWirePoints(k);
520         k++;
521
522         x2 = aPnt.X();
523         y2 = aPnt.Y();
524       }
525       // Build segment (bi-point)
526       if (j == 1)
527       {
528         aXstart = x2;
529         aYstart = y2;
530       }
531       else
532       {
533         Standard_Real *aCoordinates1 = ((Standard_Real*)(aBiPoint->Coordinates()));
534         aBiPoint++;
535
536         aCoordinates1[0] = x1;
537         aCoordinates1[1] = y1;
538         aCoordinates1[2] = x2;
539         aCoordinates1[3] = y2;
540         aCoordinates1[4] = x2 - x1;
541         aCoordinates1[5] = y2 - y1;
542       }
543       x1 = x2;
544       y1 = y2;
545     }
546   }
547
548   // Search the intersection
549   // Explore first wire
550   Standard_Integer ikEnd = 0;
551   for(Standard_Integer i = 1; i <= aNbWires; i++)
552   {
553     Standard_Integer ik = ikEnd + 1;
554     ikEnd += aWireLength(i);
555
556     // Explore second wire
557     for (Standard_Integer j = i; j <= aNbWires; j++)
558     {
559       if ( checkWiresIntersection(i, j, &ik, ikEnd, aWireLength, aBiPoints) )
560       {
561         myState = BRepMesh_SelfIntersectingWire;
562         return;
563       }
564     }
565   }
566
567   // Find holes
568   for (Standard_Integer i = aNbWires; i >= 1; i--)
569   {
570     aNbBiPoint = aWirePoints.Length() - aWireLength(i) + 1;
571     aWirePoints.Split(aNbBiPoint, aWire);
572     AnalizeWire(aWire, theUmin, theUmax, theVmin, theVmax);
573   }
574 }
575
576
577 //=======================================================================
578 //function : Perform
579 //purpose  : 
580 //=======================================================================
581 TopAbs_State BRepMesh_Classifier::Perform(const gp_Pnt2d& thePoint) const
582 {
583   Standard_Boolean isOut = Standard_False;
584   Standard_Integer aNb   = myTabClass.Length();
585   
586   for (Standard_Integer i = 1; i <= aNb; i++)
587   {
588     Standard_Integer aCur = ((CSLib_Class2d*)myTabClass(i))->SiDans(thePoint);
589     if (aCur == 0)
590     {
591       // Point is ON, but mark it as OUT
592       isOut = Standard_True;
593     }
594     else
595       isOut = myTabOrient(i)? (aCur == -1) : (aCur == 1);
596     
597     if (isOut)
598       return TopAbs_OUT;
599   }
600
601   return TopAbs_IN;
602 }
603
604
605 //=======================================================================
606 //function : Destroy
607 //purpose  : 
608 //=======================================================================
609 void BRepMesh_Classifier::Destroy()
610 {
611   Standard_Integer aNb = myTabClass.Length();
612   for (Standard_Integer i = 1; i <= aNb; i++)
613   {
614     if (myTabClass(i))
615     {
616       delete ((CSLib_Class2d*)myTabClass(i));
617       myTabClass(i) = NULL;
618     }
619   }
620 }