0026431: Can't cut a sphere from a cylinder
[occt.git] / src / IntPatch / IntPatch_WLineTool.cxx
1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
2 //
3 // This file is part of Open CASCADE Technology software library.
4 //
5 // This library is free software; you can redistribute it and/or modify it under
6 // the terms of the GNU Lesser General Public License version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
10 //
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
13
14 #include <IntPatch_WLineTool.hxx>
15
16 #include <Adaptor3d_HSurface.hxx>
17 #include <Adaptor3d_TopolTool.hxx>
18
19 //=======================================================================
20 //function : MinMax
21 //purpose  : Replaces theParMIN = MIN(theParMIN, theParMAX),
22 //                    theParMAX = MAX(theParMIN, theParMAX).
23 //
24 //           Static subfunction in IsSeamOrBound.
25 //=======================================================================
26 static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX)
27 {
28   if(theParMIN > theParMAX)
29   {
30     const Standard_Real aTmp = theParMAX;
31     theParMAX = theParMIN;
32     theParMIN = aTmp;
33   }
34 }
35
36 //=========================================================================
37 // function : FillPointsHash
38 // purpose  : Fill points hash by input data.
39 //            Static subfunction in ComputePurgedWLine.
40 //=========================================================================
41 static void FillPointsHash(const Handle(IntPatch_WLine)         &theWLine,
42                            NCollection_Array1<Standard_Integer> &thePointsHash)
43 {
44   // 1 - Delete point.
45   // 0 - Store point.
46   // -1 - Vertex point (not delete).
47   Standard_Integer i, v;
48
49   for(i = 1; i <= theWLine->NbPnts(); i++)
50     thePointsHash.SetValue(i, 0);
51
52   for(v = 1; v <= theWLine->NbVertex(); v++) 
53   {
54     IntPatch_Point aVertex = theWLine->Vertex(v);
55     Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
56     thePointsHash.SetValue(avertexindex, -1);
57   }
58 }
59
60 //=========================================================================
61 // function : MakeNewWLine
62 // purpose  : Makes new walking line according to the points hash
63 //            Static subfunction in ComputePurgedWLine and DeleteOuter.
64 //=========================================================================
65 static Handle(IntPatch_WLine) MakeNewWLine(const Handle(IntPatch_WLine)         &theWLine,
66                                            const NCollection_Array1<Standard_Integer> &thePointsHash)
67 {
68   Standard_Integer i;
69
70   Handle(IntSurf_LineOn2S) aPurgedLineOn2S = new IntSurf_LineOn2S();
71   Handle(IntPatch_WLine) aLocalWLine = new IntPatch_WLine(aPurgedLineOn2S, Standard_False);
72   Standard_Integer anOldLineIdx = 1, aVertexIdx = 1;
73   for(i = 1; i <= thePointsHash.Upper(); i++)
74   {
75     if (thePointsHash(i) == 0)
76     {
77       // Store this point.
78       aPurgedLineOn2S->Add(theWLine->Point(i));
79       anOldLineIdx++;
80     }
81     else if (thePointsHash(i) == -1)
82     {
83       // Add vertex.
84       IntPatch_Point aVertex = theWLine->Vertex(aVertexIdx++);
85       aVertex.SetParameter(anOldLineIdx++);
86       aLocalWLine->AddVertex(aVertex);
87       aPurgedLineOn2S->Add(theWLine->Point(i));
88     }
89   }
90
91   return aLocalWLine;
92 }
93
94 //=========================================================================
95 // function : MovePoint
96 // purpose  : Move point into surface param space. No interpolation used 
97 //            because walking algorithm should care for closeness to the param space.
98 //            Static subfunction in ComputePurgedWLine.
99 //=========================================================================
100 static void MovePoint(const Handle(Adaptor3d_HSurface)   &theS1,
101                       Standard_Real &U1, Standard_Real &V1)
102 {
103   if (U1 < theS1->FirstUParameter())
104     U1 = theS1->FirstUParameter();
105
106   if (U1 > theS1->LastUParameter())
107     U1 = theS1->LastUParameter();
108
109   if (V1 < theS1->FirstVParameter())
110     V1 = theS1->FirstVParameter();
111
112   if (V1 > theS1->LastVParameter())
113    V1 = theS1->LastVParameter();
114 }
115
116 //=========================================================================
117 // function : DeleteOuterPoints
118 // purpose  : Check and delete out of bounds points on walking line.
119 //            Static subfunction in ComputePurgedWLine.
120 //=========================================================================
121 static Handle(IntPatch_WLine)
122   DeleteOuterPoints(const Handle(IntPatch_WLine)       &theWLine,
123                     const Handle(Adaptor3d_HSurface)   &theS1,
124                     const Handle(Adaptor3d_HSurface)   &theS2,
125                     const Handle(Adaptor3d_TopolTool)  &theDom1,
126                     const Handle(Adaptor3d_TopolTool)  &theDom2)
127 {
128   Standard_Integer i;
129
130   NCollection_Array1<Standard_Integer> aDelOuterPointsHash(1, theWLine->NbPnts());
131   FillPointsHash(theWLine, aDelOuterPointsHash);
132
133   if (theS1->IsUPeriodic() || theS1->IsVPeriodic() ||
134       theS2->IsUPeriodic() || theS2->IsVPeriodic() )
135       return theWLine;
136
137   gp_Pnt2d aPntOnF1, aPntOnF2;
138   Standard_Real aX1, aY1, aX2, aY2;
139
140   // Iterate over points in walking line and delete which are out of bounds.
141   // Forward.
142   Standard_Boolean isAllDeleted = Standard_True;
143   Standard_Boolean aChangedFirst = Standard_False;
144   Standard_Integer aFirstGeomIdx = 1;
145   for(i = 1; i <= theWLine->NbPnts(); i++)
146   {
147     theWLine->Point(i).Parameters(aX1, aY1, aX2, aY2);
148     aPntOnF1.SetCoord(aX1, aY1);
149     aPntOnF2.SetCoord(aX2, aY2);
150
151     TopAbs_State aState1 = theDom1->Classify(aPntOnF1, Precision::Confusion());
152     TopAbs_State aState2 = theDom2->Classify(aPntOnF2, Precision::Confusion());
153
154     if (aState1 == TopAbs_OUT ||
155         aState2 == TopAbs_OUT )
156     {
157       aDelOuterPointsHash(i) = 1;
158       aChangedFirst = Standard_True;
159     }
160     else
161     {
162       isAllDeleted = Standard_False;
163
164       aFirstGeomIdx = Max (i - 1, 1);
165       if (aDelOuterPointsHash(i) == -1)
166         aFirstGeomIdx = i; // Use data what lies in (i) point / vertex.
167
168       aDelOuterPointsHash(i) = -1;
169       break;
170     }
171   }
172
173   if (isAllDeleted)
174   {
175     // ALL points are out of bounds:
176     // case boolean bcut_complex F5 and similar.
177     return theWLine;
178   }
179
180   // Backward.
181   Standard_Boolean aChangedLast = Standard_False;
182   Standard_Integer aLastGeomIdx = theWLine->NbPnts();
183   for(i = theWLine->NbPnts(); i >= 1; i--)
184   {
185     theWLine->Point(i).Parameters(aX1, aY1, aX2, aY2);
186     aPntOnF1.SetCoord(aX1, aY1);
187     aPntOnF2.SetCoord(aX2, aY2);
188
189     TopAbs_State aState1 = theDom1->Classify(aPntOnF1, Precision::Confusion());
190     TopAbs_State aState2 = theDom2->Classify(aPntOnF2, Precision::Confusion());
191
192     if (aState1 == TopAbs_OUT ||
193         aState2 == TopAbs_OUT )
194     {
195       aDelOuterPointsHash(i) = 1;
196       aChangedLast = Standard_True; // Move vertex to first good point
197     }
198     else
199     {
200       aLastGeomIdx = Min (i + 1, theWLine->NbPnts());
201       if (aDelOuterPointsHash(i) == -1)
202         aLastGeomIdx = i; // Use data what lies in (i) point / vertex.
203
204       aDelOuterPointsHash(i) = -1;
205       break;
206     }
207   }
208
209   if (!aChangedFirst && !aChangedLast)
210   {
211     // Nothing is done, return input.
212     return theWLine;
213   }
214
215   // Build new line and modify geometry of necessary vertexes.
216   Handle(IntPatch_WLine) aLocalWLine = MakeNewWLine(theWLine, aDelOuterPointsHash);
217
218   if (aChangedFirst)
219   {
220     // Vertex geometry.
221     IntPatch_Point aVertex = aLocalWLine->Vertex(1);
222     aVertex.SetValue(theWLine->Point(aFirstGeomIdx).Value());
223     Standard_Real aU1, aU2, aV1, aV2;
224     theWLine->Point(aFirstGeomIdx).Parameters(aU1, aV1, aU2, aV2);
225     MovePoint(theS1, aU1, aV1);
226     MovePoint(theS2, aU2, aV2);
227     aVertex.SetParameters(aU1, aV1, aU2, aV2);
228     aLocalWLine->Replace(1, aVertex);
229     // Change point in walking line.
230     aLocalWLine->SetPoint(1, aVertex);
231   }
232
233   if (aChangedLast)
234   {
235     // Vertex geometry.
236     IntPatch_Point aVertex = aLocalWLine->Vertex(aLocalWLine->NbVertex());
237     aVertex.SetValue(theWLine->Point(aLastGeomIdx).Value());
238     Standard_Real aU1, aU2, aV1, aV2;
239     theWLine->Point(aLastGeomIdx).Parameters(aU1, aV1, aU2, aV2);
240     MovePoint(theS1, aU1, aV1);
241     MovePoint(theS2, aU2, aV2);
242     aVertex.SetParameters(aU1, aV1, aU2, aV2);
243     aLocalWLine->Replace(aLocalWLine->NbVertex(), aVertex);
244     // Change point in walking line.
245     aLocalWLine->SetPoint(aLocalWLine->NbPnts(), aVertex);
246   }
247
248
249   return aLocalWLine;
250 }
251
252 //=========================================================================
253 // function : IsInsideIn2d
254 // purpose  : Check if aNextPnt lies inside of tube build on aBasePnt and aBaseVec.
255 //            In 2d space. Static subfunction in DeleteByTube.
256 //=========================================================================
257 static Standard_Boolean IsInsideIn2d(const gp_Pnt2d& aBasePnt,
258                                      const gp_Vec2d& aBaseVec,
259                                      const gp_Pnt2d& aNextPnt,
260                                      const Standard_Real aSquareMaxDist)
261 {
262   gp_Vec2d aVec2d(aBasePnt, aNextPnt);
263
264   //d*d = (basevec^(nextpnt-basepnt))**2 / basevec**2
265   Standard_Real aCross = aVec2d.Crossed(aBaseVec);
266   Standard_Real aSquareDist = aCross * aCross
267                             / aBaseVec.SquareMagnitude();
268
269   return (aSquareDist <= aSquareMaxDist);
270 }
271
272 //=========================================================================
273 // function : IsInsideIn3d
274 // purpose  : Check if aNextPnt lies inside of tube build on aBasePnt and aBaseVec.
275 //            In 3d space. Static subfunction in DeleteByTube.
276 //=========================================================================
277 static Standard_Boolean IsInsideIn3d(const gp_Pnt& aBasePnt,
278                                      const gp_Vec& aBaseVec,
279                                      const gp_Pnt& aNextPnt,
280                                      const Standard_Real aSquareMaxDist)
281 {
282   gp_Vec aVec(aBasePnt, aNextPnt);
283
284   //d*d = (basevec^(nextpnt-basepnt))**2 / basevec**2
285   Standard_Real aSquareDist = aVec.CrossSquareMagnitude(aBaseVec)
286                             / aBaseVec.SquareMagnitude();
287
288   return (aSquareDist <= aSquareMaxDist);
289 }
290
291 static const Standard_Integer aMinNbBadDistr = 15;
292 static const Standard_Integer aNbSingleBezier = 30;
293
294 //=========================================================================
295 // function : DeleteByTube
296 // purpose  : Check and delete points using tube criteria.
297 //            Static subfunction in ComputePurgedWLine.
298 //=========================================================================
299 static Handle(IntPatch_WLine)
300   DeleteByTube(const Handle(IntPatch_WLine)       &theWLine,
301                const Handle(Adaptor3d_HSurface)   &theS1,
302                const Handle(Adaptor3d_HSurface)   &theS2)
303 {
304   // III: Check points for tube criteria:
305   // Workaround to handle case of small amount points after purge.
306   // Test "boolean boptuc_complex B5" and similar.
307   Standard_Integer aNbPnt = 0 , i;
308
309   if (theWLine->NbPnts() <= 2)
310     return theWLine;
311
312   NCollection_Array1<Standard_Integer> aNewPointsHash(1, theWLine->NbPnts());
313   FillPointsHash(theWLine, aNewPointsHash);
314   
315   // Inital computations.
316   Standard_Real UonS1[3], VonS1[3], UonS2[3], VonS2[3];
317   theWLine->Point(1).ParametersOnS1(UonS1[0], VonS1[0]);
318   theWLine->Point(2).ParametersOnS1(UonS1[1], VonS1[1]);
319   theWLine->Point(1).ParametersOnS2(UonS2[0], VonS2[0]);
320   theWLine->Point(2).ParametersOnS2(UonS2[1], VonS2[1]);
321
322   gp_Pnt2d aBase2dPnt1(UonS1[0], VonS1[0]);
323   gp_Pnt2d aBase2dPnt2(UonS2[0], VonS2[0]);
324   gp_Vec2d aBase2dVec1(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0]);
325   gp_Vec2d aBase2dVec2(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0]);
326   gp_Pnt   aBase3dPnt = theWLine->Point(1).Value();
327   gp_Vec   aBase3dVec(theWLine->Point(1).Value(), theWLine->Point(2).Value());
328
329   // Choose base tolerance and scale it to pipe algorithm.
330   const Standard_Real aBaseTolerance = Precision::Approximation();
331   Standard_Real aResS1Tol = Min(theS1->UResolution(aBaseTolerance),
332                                 theS1->VResolution(aBaseTolerance));
333   Standard_Real aResS2Tol = Min(theS2->UResolution(aBaseTolerance),
334                                 theS2->VResolution(aBaseTolerance));
335   Standard_Real aTol1 = aResS1Tol * aResS1Tol;
336   Standard_Real aTol2 = aResS2Tol * aResS2Tol;
337   Standard_Real aTol3d = aBaseTolerance * aBaseTolerance;
338
339   const Standard_Real aLimitCoeff = 0.99 * 0.99;
340   for(i = 3; i <= theWLine->NbPnts(); i++)
341   {
342     Standard_Boolean isDeleteState = Standard_False;
343
344     theWLine->Point(i).ParametersOnS1(UonS1[2], VonS1[2]);
345     theWLine->Point(i).ParametersOnS2(UonS2[2], VonS2[2]);
346     gp_Pnt2d aPnt2dOnS1(UonS1[2], VonS1[2]);
347     gp_Pnt2d aPnt2dOnS2(UonS2[2], VonS2[2]);
348     const gp_Pnt& aPnt3d = theWLine->Point(i).Value();
349
350     if (aNewPointsHash(i - 1) != - 1 &&
351         IsInsideIn2d(aBase2dPnt1, aBase2dVec1, aPnt2dOnS1, aTol1) &&
352         IsInsideIn2d(aBase2dPnt2, aBase2dVec2, aPnt2dOnS2, aTol2) &&
353         IsInsideIn3d(aBase3dPnt, aBase3dVec, aPnt3d, aTol3d) )
354     {
355       // Handle possible uneven parametrization on one of 2d subspaces.
356       // Delete point only when expected lengths are close to each other (aLimitCoeff).
357       // Example:
358       // c2d1 - line
359       // c3d - line
360       // c2d2 - geometrically line, but have uneven parametrization -> c2d2 is bspline.
361       gp_XY aPntOnS1[2]= { gp_XY(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0])
362                          , gp_XY(UonS1[2] - UonS1[1], VonS1[2] - VonS1[1])};
363       gp_XY aPntOnS2[2]= { gp_XY(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0])
364                          , gp_XY(UonS2[2] - UonS2[1], VonS2[2] - VonS2[1])};
365
366       Standard_Real aStepOnS1 = aPntOnS1[0].SquareModulus() / aPntOnS1[1].SquareModulus();
367       Standard_Real aStepOnS2 = aPntOnS2[0].SquareModulus() / aPntOnS2[1].SquareModulus();
368
369       Standard_Real aStepCoeff = Min(aStepOnS1, aStepOnS2) / Max(aStepOnS1, aStepOnS2);
370
371       if (aStepCoeff > aLimitCoeff)
372       {
373         // Set hash flag to "Delete" state.
374         isDeleteState = Standard_True;
375         aNewPointsHash.SetValue(i - 1, 1);
376
377         // Change middle point.
378         UonS1[1] = UonS1[2];
379         UonS2[1] = UonS2[2];
380         VonS1[1] = VonS1[2];
381         VonS2[1] = VonS2[2];
382       }
383     }
384
385     if (!isDeleteState)
386     {
387       // Compute new pipe parameters.
388       UonS1[0] = UonS1[1];
389       VonS1[0] = VonS1[1];
390       UonS2[0] = UonS2[1];
391       VonS2[0] = VonS2[1];
392
393       UonS1[1] = UonS1[2];
394       VonS1[1] = VonS1[2];
395       UonS2[1] = UonS2[2];
396       VonS2[1] = VonS2[2];
397
398       aBase2dPnt1.SetCoord(UonS1[0], VonS1[0]);
399       aBase2dPnt2.SetCoord(UonS2[0], VonS2[0]);
400       aBase2dVec1.SetCoord(UonS1[1] - UonS1[0], VonS1[1] - VonS1[0]);
401       aBase2dVec2.SetCoord(UonS2[1] - UonS2[0], VonS2[1] - VonS2[0]);
402       aBase3dPnt = theWLine->Point(i - 1).Value();
403       aBase3dVec = gp_Vec(theWLine->Point(i - 1).Value(), theWLine->Point(i).Value());
404
405       aNbPnt++;
406     }
407   }
408
409   // Workaround to handle case of small amount of points after purge.
410   // Test "boolean boptuc_complex B5" and similar.
411   // This is possible since there are at least two points.
412   if (aNewPointsHash(1) == -1 &&
413       aNewPointsHash(2) == -1 &&
414       aNbPnt <= 3)
415   {
416     // Delete first.
417     aNewPointsHash(1) = 1;
418   }
419   if (aNewPointsHash(theWLine->NbPnts() - 1) == -1 &&
420       aNewPointsHash(theWLine->NbPnts()    ) == -1 &&
421       aNbPnt <= 3)
422   {
423     // Delete last.
424     aNewPointsHash(theWLine->NbPnts()) = 1;
425   }
426
427   // Purgre when too small amount of points left.
428   if (aNbPnt <= 2)
429   {
430     for(i = aNewPointsHash.Lower(); i <= aNewPointsHash.Upper(); i++)
431     {
432       if (aNewPointsHash(i) != -1)
433       {
434         aNewPointsHash(i) = 1;
435       }
436     }
437   }
438
439   // Handle possible bad distribution of points, 
440   // which are will converted into one single bezier curve (less than 30 points).
441   // Make distribution more even:
442   // max step will be nearly to 0.1 of param distance.
443   if (aNbPnt + 2 > aMinNbBadDistr &&
444       aNbPnt + 2 < aNbSingleBezier )
445   {
446     for(Standard_Integer anIdx = 1; anIdx <= 8; anIdx++)
447     {
448       Standard_Integer aHashIdx = 
449         Standard_Integer(anIdx * theWLine->NbPnts() / 9);
450
451       //Store this point.
452       aNewPointsHash(aHashIdx) = 0;
453     }
454   }
455
456   return MakeNewWLine(theWLine, aNewPointsHash);
457 }
458
459 //=======================================================================
460 //function : IsSeam
461 //purpose  : Returns:
462 //            0 - if interval [theU1, theU2] does not intersect the "seam-edge"
463 //                or if "seam-edge" do not exist;
464 //            1 - if interval (theU1, theU2) intersect the "seam-edge".
465 //            2 - if theU1 or/and theU2 lie ON the "seam-edge"
466 //
467 //ATTENTION!!!
468 //  If (theU1 == theU2) then this function will return only both 0 or 2.
469 //
470 //           Static subfunction in IsSeamOrBound.
471 //=======================================================================
472 static Standard_Integer IsSeam( const Standard_Real theU1,
473                                 const Standard_Real theU2,
474                                 const Standard_Real thePeriod)
475 {
476   if(IsEqual(thePeriod, 0.0))
477     return 0;
478
479   //If interval [theU1, theU2] intersect seam-edge then there exists an integer
480   //number N such as 
481   //    (theU1 <= T*N <= theU2) <=> (theU1/T <= N <= theU2/T),
482   //where T is the period.
483   //I.e. the inerval [theU1/T, theU2/T] must contain at least one
484   //integer number. In this case, Floor(theU1/T) and Floor(theU2/T)
485   //return different values or theU1/T is strictly integer number.
486   //Examples:
487   //  1. theU1/T==2.8, theU2/T==3.5 => Floor(theU1/T) == 2, Floor(theU2/T) == 3.
488   //  2. theU1/T==2.0, theU2/T==2.6 => Floor(theU1/T) == Floor(theU2/T) == 2.
489
490   const Standard_Real aVal1 = theU1/thePeriod,
491                       aVal2 = theU2/thePeriod;
492   const Standard_Integer aPar1 = static_cast<Standard_Integer>(Floor(aVal1));
493   const Standard_Integer aPar2 = static_cast<Standard_Integer>(Floor(aVal2));
494   if(aPar1 != aPar2)
495   {//Interval (theU1, theU2] intersects seam-edge
496     if(IsEqual(aVal2, static_cast<Standard_Real>(aPar2)))
497     {//aVal2 is an integer number => theU2 lies ON the "seam-edge"
498       return 2;
499     }
500
501     return 1;
502   }
503
504   //Here, aPar1 == aPar2. 
505
506   if(IsEqual(aVal1, static_cast<Standard_Real>(aPar1)))
507   {//aVal1 is an integer number => theU1 lies ON the "seam-edge"
508     return 2;
509   }
510
511   //If aVal2 is a true integer number then always (aPar1 != aPar2).
512
513   return 0;
514 }
515
516 //=======================================================================
517 //function : IsSeamOrBound
518 //purpose  : Returns TRUE if segment [thePtf, thePtl] intersects "seam-edge"
519 //            (if it exist) or surface boundaries and both thePtf and thePtl do
520 //            not match "seam-edge" or boundaries.
521 //           Point thePtmid lies in this segment. If thePtmid match
522 //            "seam-edge" or boundaries strictly (without any tolerance) then
523 //            the function will return TRUE.
524 //            See comments in function body for detail information.
525 //
526 //           Static subfunction in JoinWLines.
527 //=======================================================================
528 static Standard_Boolean IsSeamOrBound(const IntSurf_PntOn2S& thePtf,
529                                       const IntSurf_PntOn2S& thePtl,
530                                       const IntSurf_PntOn2S& thePtmid,
531                                       const Standard_Real theU1Period,
532                                       const Standard_Real theU2Period,
533                                       const Standard_Real theV1Period,
534                                       const Standard_Real theV2Period,
535                                       const Standard_Real theUfSurf1,
536                                       const Standard_Real theUlSurf1,
537                                       const Standard_Real theVfSurf1,
538                                       const Standard_Real theVlSurf1,
539                                       const Standard_Real theUfSurf2,
540                                       const Standard_Real theUlSurf2,
541                                       const Standard_Real theVfSurf2,
542                                       const Standard_Real theVlSurf2)
543 {
544   Standard_Real aU11 = 0.0, aU12 = 0.0, aV11 = 0.0, aV12 = 0.0;
545   Standard_Real aU21 = 0.0, aU22 = 0.0, aV21 = 0.0, aV22 = 0.0;
546   thePtf.Parameters(aU11, aV11, aU12, aV12);
547   thePtl.Parameters(aU21, aV21, aU22, aV22);
548
549   MinMax(aU11, aU21);
550   MinMax(aV11, aV21);
551   MinMax(aU12, aU22);
552   MinMax(aV12, aV22);
553
554   if((aU11 - theUfSurf1)*(aU21 - theUfSurf1) < 0.0)
555   {//Interval [aU11, aU21] intersects theUfSurf1
556     return Standard_True;
557   }
558
559   if((aU11 - theUlSurf1)*(aU21 - theUlSurf1) < 0.0)
560   {//Interval [aU11, aU21] intersects theUlSurf1
561     return Standard_True;
562   }
563
564   if((aV11 - theVfSurf1)*(aV21 - theVfSurf1) < 0.0)
565   {//Interval [aV11, aV21] intersects theVfSurf1
566     return Standard_True;
567   }
568
569   if((aV11 - theVlSurf1)*(aV21 - theVlSurf1) < 0.0)
570   {//Interval [aV11, aV21] intersects theVlSurf1
571     return Standard_True;
572   }
573
574   if((aU12 - theUfSurf2)*(aU22 - theUfSurf2) < 0.0)
575   {//Interval [aU12, aU22] intersects theUfSurf2
576     return Standard_True;
577   }
578
579   if((aU12 - theUlSurf2)*(aU22 - theUlSurf2) < 0.0)
580   {//Interval [aU12, aU22] intersects theUlSurf2
581     return Standard_True;
582   }
583
584   if((aV12 - theVfSurf2)*(aV22 - theVfSurf2) < 0.0)
585   {//Interval [aV12, aV22] intersects theVfSurf2
586     return Standard_True;
587   }
588
589   if((aV12 - theVlSurf2)*(aV22 - theVlSurf2) < 0.0)
590   {//Interval [aV12, aV22] intersects theVlSurf2
591     return Standard_True;
592   }
593
594   if(IsSeam(aU11, aU21, theU1Period))
595     return Standard_True;
596
597   if(IsSeam(aV11, aV21, theV1Period))
598     return Standard_True;
599
600   if(IsSeam(aU12, aU22, theU2Period))
601     return Standard_True;
602
603   if(IsSeam(aV12, aV22, theV2Period))
604     return Standard_True;
605
606   /*
607     The segment [thePtf, thePtl] does not intersect the boundaries and
608     the seam-edge of the surfaces.
609     Nevertheless, following situation is possible:
610
611                   seam or
612                    bound
613                      |
614         thePtf  *    |
615                      |
616                      * thePtmid
617           thePtl  *  |
618                      |
619
620     This case must be processed, too.
621   */
622
623   Standard_Real aU1 = 0.0, aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
624   thePtmid.Parameters(aU1, aV1, aU2, aV2);
625
626   if(IsEqual(aU1, theUfSurf1) || IsEqual(aU1, theUlSurf1))
627     return Standard_True;
628
629   if(IsEqual(aU2, theUfSurf2) || IsEqual(aU2, theUlSurf2))
630     return Standard_True;
631
632   if(IsEqual(aV1, theVfSurf1) || IsEqual(aV1, theVlSurf1))
633     return Standard_True;
634
635   if(IsEqual(aV2, theVfSurf2) || IsEqual(aV2, theVlSurf2))
636     return Standard_True;
637
638   if(IsSeam(aU1, aU1, theU1Period))
639     return Standard_True;
640
641   if(IsSeam(aU2, aU2, theU2Period))
642     return Standard_True;
643
644   if(IsSeam(aV1, aV1, theV1Period))
645     return Standard_True;
646
647   if(IsSeam(aV2, aV2, theV2Period))
648     return Standard_True;
649
650   return Standard_False;
651 }
652
653 //=========================================================================
654 // function : ComputePurgedWLine
655 // purpose  :
656 //=========================================================================
657 Handle(IntPatch_WLine) IntPatch_WLineTool::
658   ComputePurgedWLine(const Handle(IntPatch_WLine)       &theWLine,
659                      const Handle(Adaptor3d_HSurface)   &theS1,
660                      const Handle(Adaptor3d_HSurface)   &theS2,
661                      const Handle(Adaptor3d_TopolTool)  &theDom1,
662                      const Handle(Adaptor3d_TopolTool)  &theDom2)
663 {
664   Standard_Integer i, k, v, nb, nbvtx;
665   Handle(IntPatch_WLine) aResult;
666   nbvtx = theWLine->NbVertex();
667   nb = theWLine->NbPnts();
668   if (nb==2)
669   {
670     const IntSurf_PntOn2S& p1 = theWLine->Point(1);
671     const IntSurf_PntOn2S& p2 = theWLine->Point(2);
672     if(p1.Value().IsEqual(p2.Value(), gp::Resolution()))
673       return aResult;
674   }
675
676   Handle(IntPatch_WLine) aLocalWLine;
677   Handle(IntPatch_WLine) aTmpWLine = theWLine;
678   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
679   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
680   for(i = 1; i <= nb; i++)
681     aLineOn2S->Add(theWLine->Point(i));
682
683   for(v = 1; v <= nbvtx; v++)
684     aLocalWLine->AddVertex(theWLine->Vertex(v));
685
686   // I: Delete equal points
687   for(i = 1; i <= aLineOn2S->NbPoints(); i++)
688   {
689     Standard_Integer aStartIndex = i + 1;
690     Standard_Integer anEndIndex = i + 5;
691     nb = aLineOn2S->NbPoints();
692     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
693
694     if((aStartIndex > nb) || (anEndIndex <= 1))
695       continue;
696
697     k = aStartIndex;
698
699     while(k <= anEndIndex)
700     {
701       if(i != k)
702       {
703         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
704         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
705         
706         Standard_Real UV[8];
707         p1.Parameters(UV[0], UV[1], UV[2], UV[3]);
708         p2.Parameters(UV[4], UV[5], UV[6], UV[7]);
709
710         Standard_Real aMax = Abs(UV[0]);
711         for(Standard_Integer anIdx = 1; anIdx < 8; anIdx++)
712         {
713           if (aMax < Abs(UV[anIdx]))
714             aMax = Abs(UV[anIdx]);
715         }
716
717         if(p1.Value().IsEqual(p2.Value(), gp::Resolution()) ||
718            Abs(UV[0] - UV[4]) + Abs(UV[1] - UV[5]) < 1.0e-16 * aMax ||
719            Abs(UV[2] - UV[6]) + Abs(UV[3] - UV[7]) < 1.0e-16 * aMax )
720         {
721           aTmpWLine = aLocalWLine;
722           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
723           
724           for(v = 1; v <= aTmpWLine->NbVertex(); v++)
725           {
726             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
727             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
728
729             if(avertexindex >= k)
730             {
731               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
732             }
733             aLocalWLine->AddVertex(aVertex);
734           }
735           aLineOn2S->RemovePoint(k);
736           anEndIndex--;
737           continue;
738         }
739       }
740       k++;
741     }
742   }
743
744   if (aLineOn2S->NbPoints() <= 2)
745   {
746     if (aLineOn2S->NbPoints() == 2)
747       return aLocalWLine;
748     else
749       return aResult;
750   }
751
752   // Avoid purge in case of C0 continuity:
753   // Intersection approximator may produce invalid curve after purge, example:
754   // bugs modalg_5 bug24731.
755   // Do not run purger when base number of points is too small.
756   if (theS1->UContinuity() == GeomAbs_C0 ||
757       theS1->VContinuity() == GeomAbs_C0 ||
758       theS2->UContinuity() == GeomAbs_C0 ||
759       theS2->VContinuity() == GeomAbs_C0 ||
760       nb < aNbSingleBezier)
761   {
762     return aLocalWLine;
763   }
764
765   // II: Delete out of borders points.
766   Handle(IntPatch_WLine) aLocalWLineOuter = 
767     DeleteOuterPoints(aLocalWLine, theS1, theS2, theDom1, theDom2);
768
769   // III: Delete points by tube criteria.
770   Handle(IntPatch_WLine) aLocalWLineTube = 
771     DeleteByTube(aLocalWLineOuter, theS1, theS2);
772
773   if(aLocalWLineTube->NbPnts() > 1)
774   {
775     aResult = aLocalWLineTube;
776   }
777   return aResult;
778 }
779
780
781 //=======================================================================
782 //function : JoinWLines
783 //purpose  :
784 //=======================================================================
785 void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin,
786                                     IntPatch_SequenceOfPoint& theSPnt,
787                                     const Standard_Real theTol3D,
788                                     const Standard_Real theU1Period,
789                                     const Standard_Real theU2Period,
790                                     const Standard_Real theV1Period,
791                                     const Standard_Real theV2Period,
792                                     const Standard_Real theUfSurf1,
793                                     const Standard_Real theUlSurf1,
794                                     const Standard_Real theVfSurf1,
795                                     const Standard_Real theVlSurf1,
796                                     const Standard_Real theUfSurf2,
797                                     const Standard_Real theUlSurf2,
798                                     const Standard_Real theVfSurf2,
799                                     const Standard_Real theVlSurf2)
800 {
801   if(theSlin.Length() == 0)
802     return;
803
804   for(Standard_Integer aNumOfLine1 = 1; aNumOfLine1 <= theSlin.Length(); aNumOfLine1++)
805   {
806     Handle(IntPatch_WLine) aWLine1 (Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNumOfLine1)));
807
808     if(aWLine1.IsNull())
809     {//We must have failed to join not-point-lines
810       continue;
811     }
812
813     const Standard_Integer aNbPntsWL1 = aWLine1->NbPnts();
814     const IntSurf_PntOn2S& aPntFW1 = aWLine1->Point(1);
815     const IntSurf_PntOn2S& aPntLW1 = aWLine1->Point(aNbPntsWL1);
816
817     for(Standard_Integer aNPt = 1; aNPt <= theSPnt.Length(); aNPt++)
818     {
819       const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNPt).PntOn2S();
820
821       if( aPntCur.IsSame(aPntFW1, Precision::Confusion()) ||
822         aPntCur.IsSame(aPntLW1, Precision::Confusion()))
823       {
824         theSPnt.Remove(aNPt);
825         aNPt--;
826       }
827     }
828
829     Standard_Boolean hasBeenRemoved = Standard_False;
830     for(Standard_Integer aNumOfLine2 = aNumOfLine1 + 1; aNumOfLine2 <= theSlin.Length(); aNumOfLine2++)
831     {
832       Handle(IntPatch_WLine) aWLine2 (Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNumOfLine2)));
833
834       if(aWLine2.IsNull())
835         continue;
836
837       const Standard_Integer aNbPntsWL2 = aWLine2->NbPnts();
838
839       const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
840       const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aNbPntsWL1);
841
842       const IntSurf_PntOn2S& aPntFWL2 = aWLine2->Point(1);
843       const IntSurf_PntOn2S& aPntLWL2 = aWLine2->Point(aNbPntsWL2);
844
845       if(aPntFWL1.IsSame(aPntFWL2, Precision::Confusion()))
846       {
847         const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2);
848         const IntSurf_PntOn2S& aPt2 = aWLine2->Point(2);
849         if(!IsSeamOrBound(aPt1, aPt2, aPntFWL1, theU1Period, theU2Period,
850                           theV1Period, theV2Period, theUfSurf1, theUlSurf1,
851                           theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
852                           theVfSurf2, theVlSurf2))
853         {
854           aWLine1->ClearVertexes();
855           for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++)
856           {
857             const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt);
858             aWLine1->Curve()->InsertBefore(1, aPt);
859           }
860
861           aWLine1->ComputeVertexParameters(theTol3D);
862
863           theSlin.Remove(aNumOfLine2);
864           aNumOfLine2--;
865           hasBeenRemoved = Standard_True;
866
867           continue;
868         }
869       }
870
871       if(aPntFWL1.IsSame(aPntLWL2, Precision::Confusion()))
872       {
873         const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2);
874         const IntSurf_PntOn2S& aPt2 = aWLine2->Point(aNbPntsWL2-1);
875         if(!IsSeamOrBound(aPt1, aPt2, aPntFWL1, theU1Period, theU2Period,
876                           theV1Period, theV2Period, theUfSurf1, theUlSurf1,
877                           theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
878                           theVfSurf2, theVlSurf2))
879         {
880           aWLine1->ClearVertexes();
881           for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--)
882           {
883             const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt);
884             aWLine1->Curve()->InsertBefore(1, aPt);
885           }
886
887           aWLine1->ComputeVertexParameters(theTol3D);
888
889           theSlin.Remove(aNumOfLine2);
890           aNumOfLine2--;
891           hasBeenRemoved = Standard_True;
892
893           continue;
894         }
895       }
896
897       if(aPntLWL1.IsSame(aPntFWL2, Precision::Confusion()))
898       {
899         const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1-1);
900         const IntSurf_PntOn2S& aPt2 = aWLine2->Point(2);
901         if(!IsSeamOrBound(aPt1, aPt2, aPntLWL1, theU1Period, theU2Period,
902                           theV1Period, theV2Period, theUfSurf1, theUlSurf1,
903                           theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
904                           theVfSurf2, theVlSurf2))
905         {
906           aWLine1->ClearVertexes();
907           for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++)
908           {
909             const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt);
910             aWLine1->Curve()->Add(aPt);
911           }
912
913           aWLine1->ComputeVertexParameters(theTol3D);
914
915           theSlin.Remove(aNumOfLine2);
916           aNumOfLine2--;
917           hasBeenRemoved = Standard_True;
918
919           continue;
920         }
921       }
922
923       if(aPntLWL1.IsSame(aPntLWL2, Precision::Confusion()))
924       {
925         const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1-1);
926         const IntSurf_PntOn2S& aPt2 = aWLine2->Point(aNbPntsWL2-1);
927         if(!IsSeamOrBound(aPt1, aPt2, aPntLWL1, theU1Period, theU2Period,
928                           theV1Period, theV2Period, theUfSurf1, theUlSurf1,
929                           theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
930                           theVfSurf2, theVlSurf2))
931         {
932           aWLine1->ClearVertexes();
933           for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--)
934           {
935             const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt);
936             aWLine1->Curve()->Add(aPt);
937           }
938
939           aWLine1->ComputeVertexParameters(theTol3D);
940
941           theSlin.Remove(aNumOfLine2);
942           aNumOfLine2--;
943           hasBeenRemoved = Standard_True;
944
945           continue;
946         }
947       }
948     }
949
950     if(hasBeenRemoved)
951       aNumOfLine1--;
952   }
953 }