0030913: Invalid result of Fusing slices
[occt.git] / src / IntTools / IntTools_EdgeEdge.cxx
1 // Created by: Eugeny MALTCHIKOV
2 // Copyright (c) 2013-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15
16 #include <Bnd_Box.hxx>
17 #include <BndLib_Add3dCurve.hxx>
18 #include <TColStd_MapOfInteger.hxx>
19 #include <BRep_Tool.hxx>
20 #include <BRepAdaptor_Curve.hxx>
21 #include <ElCLib.hxx>
22 #include <Geom_BezierCurve.hxx>
23 #include <Geom_BSplineCurve.hxx>
24 #include <Geom_Line.hxx>
25 #include <Geom_Circle.hxx>
26 #include <Geom_Curve.hxx>
27 #include <Geom_Ellipse.hxx>
28 #include <Geom_OffsetCurve.hxx>
29 #include <GeomAPI_ProjectPointOnCurve.hxx>
30 #include <gp_Dir.hxx>
31 #include <gp_Lin.hxx>
32 #include <IntTools_CommonPrt.hxx>
33 #include <IntTools_EdgeEdge.hxx>
34 #include <IntTools_Range.hxx>
35 #include <IntTools_Tools.hxx>
36 #include <TopoDS_Edge.hxx>
37 #include <TopoDS_Iterator.hxx>
38 #include <BRepExtrema_DistShapeShape.hxx>
39
40 static 
41   void BndBuildBox(const BRepAdaptor_Curve& theBAC,
42                    const Standard_Real aT1,
43                    const Standard_Real aT2,
44                    const Standard_Real theTol,
45                    Bnd_Box& theBox);
46 static
47   Standard_Real PointBoxDistance(const Bnd_Box& aB,
48                                  const gp_Pnt& aP);
49 static 
50   Standard_Integer SplitRangeOnSegments(const Standard_Real aT1, 
51                                         const Standard_Real aT2,
52                                         const Standard_Real theResolution,
53                                         const Standard_Integer theNbSeg,
54                                         IntTools_SequenceOfRanges& theSegments);
55 static
56  Standard_Integer DistPC(const Standard_Real aT1, 
57                          const Handle(Geom_Curve)& theC1,
58                          const Standard_Real theCriteria, 
59                          GeomAPI_ProjectPointOnCurve& theProjector,
60                          Standard_Real& aD, 
61                          Standard_Real& aT2,
62                          const Standard_Integer iC = 1);
63 static
64  Standard_Integer DistPC(const Standard_Real aT1, 
65                          const Handle(Geom_Curve)& theC1,
66                          const Standard_Real theCriteria,
67                          GeomAPI_ProjectPointOnCurve& theProjector, 
68                          Standard_Real& aD, 
69                          Standard_Real& aT2,
70                          Standard_Real& aDmax,
71                          Standard_Real& aT1max,
72                          Standard_Real& aT2max,
73                          const Standard_Integer iC = 1);
74 static
75   Standard_Integer FindDistPC(const Standard_Real aT1A, 
76                               const Standard_Real aT1B,
77                               const Handle(Geom_Curve)& theC1,
78                               const Standard_Real theCriteria,
79                               const Standard_Real theEps,
80                               GeomAPI_ProjectPointOnCurve& theProjector,
81                               Standard_Real& aDmax, 
82                               Standard_Real& aT1max,
83                               Standard_Real& aT2max,
84                               const Standard_Boolean bMaxDist = Standard_True);
85 static
86   Standard_Real ResolutionCoeff(const BRepAdaptor_Curve& theBAC,
87                                 const IntTools_Range& theRange);
88 static
89   Standard_Real Resolution(const Handle(Geom_Curve)& theCurve,
90                            const GeomAbs_CurveType theCurveType,
91                            const Standard_Real theResCoeff,
92                            const Standard_Real theR3D);
93 static
94   Standard_Real CurveDeflection(const BRepAdaptor_Curve& theBAC,
95                                 const IntTools_Range& theRange);
96 static 
97   Standard_Boolean IsClosed(const Handle(Geom_Curve)& theCurve,
98                             const Standard_Real aT1,
99                             const Standard_Real aT2,
100                             const Standard_Real theTol,
101                             const Standard_Real theRes);
102 static 
103   Standard_Integer TypeToInteger(const GeomAbs_CurveType theCType);
104
105 //=======================================================================
106 //function : Prepare
107 //purpose  : 
108 //=======================================================================
109 void IntTools_EdgeEdge::Prepare()
110 {
111   GeomAbs_CurveType aCT1, aCT2;
112   Standard_Integer iCT1, iCT2;
113   //
114   myCurve1.Initialize(myEdge1);
115   myCurve2.Initialize(myEdge2);
116   //
117   if (myRange1.First() == 0. && myRange1.Last() == 0.) {
118     myRange1.SetFirst(myCurve1.FirstParameter());
119     myRange1.SetLast (myCurve1.LastParameter());
120   }
121   //
122   if (myRange2.First() == 0. && myRange2.Last() == 0.) {
123     myRange2.SetFirst(myCurve2.FirstParameter());
124     myRange2.SetLast (myCurve2.LastParameter());
125   }
126   //
127   aCT1 = myCurve1.GetType();
128   aCT2 = myCurve2.GetType();
129   //
130   iCT1 = TypeToInteger(aCT1);
131   iCT2 = TypeToInteger(aCT2);
132   //
133   if (iCT1 == iCT2) {
134     if (iCT1 != 0) {
135       //compute deflection
136       Standard_Real aC1, aC2;
137       //
138       aC2 = CurveDeflection(myCurve2, myRange2);
139       aC1 = (aC2 > Precision::Confusion()) ? 
140         CurveDeflection(myCurve1, myRange1) : 1.;
141       //
142       if (aC1 < aC2) {
143         --iCT1;
144       }
145     }
146   }
147   //
148   if (iCT1 < iCT2) {
149     TopoDS_Edge tmpE = myEdge1;
150     myEdge1 = myEdge2;
151     myEdge2 = tmpE;
152     //
153     BRepAdaptor_Curve tmpC = myCurve1;
154     myCurve1 = myCurve2;
155     myCurve2 = tmpC;
156     //
157     IntTools_Range tmpR = myRange1;
158     myRange1 = myRange2;
159     myRange2 = tmpR;
160     //
161     mySwap = Standard_True;
162   }
163   //
164   Standard_Real aTolAdd = myFuzzyValue / 2.;
165   myTol1 = myCurve1.Tolerance() + aTolAdd;
166   myTol2 = myCurve2.Tolerance() + aTolAdd;
167   myTol = myTol1 + myTol2;
168   //
169   if (iCT1 != 0 || iCT2 != 0) {
170     Standard_Real f, l, aTM;
171     //
172     myGeom1 = BRep_Tool::Curve(myEdge1, f, l);
173     myGeom2 = BRep_Tool::Curve(myEdge2, f, l);
174     //
175     myResCoeff1 = ResolutionCoeff(myCurve1, myRange1);
176     myResCoeff2 = ResolutionCoeff(myCurve2, myRange2);
177     //
178     myRes1 = Resolution(myCurve1.Curve().Curve(), myCurve1.GetType(), myResCoeff1, myTol1);
179     myRes2 = Resolution(myCurve2.Curve().Curve(), myCurve2.GetType(), myResCoeff2, myTol2);
180     //
181     myPTol1 = 5.e-13;
182     aTM = Max(fabs(myRange1.First()), fabs(myRange1.Last()));
183     if (aTM > 999.) {
184       myPTol1 = 5.e-16 * aTM;
185     }
186     //
187     myPTol2 = 5.e-13;
188     aTM = Max(fabs(myRange2.First()), fabs(myRange2.Last()));
189     if (aTM > 999.) {
190       myPTol2 = 5.e-16 * aTM;
191     }
192   }
193 }
194
195 //=======================================================================
196 //function : Perform
197 //purpose  : 
198 //=======================================================================
199 void IntTools_EdgeEdge::Perform()
200 {
201   //1. Check data
202   CheckData();
203   if (myErrorStatus) {
204     return;
205   }
206   //
207   //2. Prepare Data
208   Prepare();
209   //
210   //3.1. Check Line/Line case
211   if (myCurve1.GetType() == GeomAbs_Line &&
212       myCurve2.GetType() == GeomAbs_Line) {
213     ComputeLineLine();
214     return;
215   }
216   //
217   if (myQuickCoincidenceCheck) {
218     if (IsCoincident()) {
219       Standard_Real aT11, aT12, aT21, aT22;
220       //
221       myRange1.Range(aT11, aT12);
222       myRange2.Range(aT21, aT22);
223       AddSolution(aT11, aT12, aT21, aT22, TopAbs_EDGE);
224       return;
225     }
226   }
227   //
228   if ((myCurve1.GetType() <= GeomAbs_Parabola && myCurve2.GetType() <= GeomAbs_Parabola) &&
229       (myCurve1.GetType() == GeomAbs_Line || myCurve2.GetType() == GeomAbs_Line))
230   {
231     //Improvement of performance for cases of searching common parts between line  
232     //and analytical curve. This code allows to define that edges have no
233     //common parts more fast, then regular algorithm (FindSolution(...))
234     //Check minimal distance between edges
235     BRepExtrema_DistShapeShape  aMinDist(myEdge1, myEdge2, Extrema_ExtFlag_MIN);
236     if (aMinDist.IsDone())
237     {
238       Standard_Real d = aMinDist.Value();
239       if (d > 1.1 * myTol)
240       {
241         return;
242       }
243     }
244   }
245
246   IntTools_SequenceOfRanges aRanges1, aRanges2;
247   //
248   //3.2. Find ranges containig solutions
249   Standard_Boolean bSplit2;
250   FindSolutions(aRanges1, aRanges2, bSplit2);
251   //
252   //4. Merge solutions and save common parts
253   MergeSolutions(aRanges1, aRanges2, bSplit2);
254 }
255
256 //=======================================================================
257 //function :  IsCoincident
258 //purpose  : 
259 //=======================================================================
260 Standard_Boolean IntTools_EdgeEdge::IsCoincident() 
261 {
262   Standard_Integer i, iCnt, aNbSeg, aNbP2;
263   Standard_Real dT, aT1, aCoeff, aTresh, aD;
264   Standard_Real aT11, aT12, aT21, aT22;
265   GeomAPI_ProjectPointOnCurve aProjPC;
266   gp_Pnt aP1;
267   //
268   aTresh=0.5;
269   aNbSeg=23;
270   myRange1.Range(aT11, aT12);
271   myRange2.Range(aT21, aT22);
272   //
273   aProjPC.Init(myGeom2, aT21, aT22);
274   //
275   dT=(aT12-aT11)/aNbSeg;
276   //
277   iCnt=0;
278   for(i=0; i <= aNbSeg; ++i) {
279     aT1 = aT11+i*dT;
280     myGeom1->D0(aT1, aP1);
281     //
282     aProjPC.Perform(aP1);
283     aNbP2=aProjPC.NbPoints();
284     if (!aNbP2) {
285       continue;
286     }
287     //
288     aD=aProjPC.LowerDistance();
289     if(aD < myTol) {
290       ++iCnt; 
291     }
292   }
293   //
294   aCoeff=(Standard_Real)iCnt/((Standard_Real)aNbSeg+1);
295   return aCoeff > aTresh;
296 }
297 //=======================================================================
298 //function : FindSolutions
299 //purpose  : 
300 //=======================================================================
301 void IntTools_EdgeEdge::FindSolutions(IntTools_SequenceOfRanges& theRanges1,
302                                       IntTools_SequenceOfRanges& theRanges2,
303                                       Standard_Boolean& bSplit2)
304 {
305   Standard_Boolean bIsClosed2;
306   Standard_Real aT11, aT12, aT21, aT22;
307   Bnd_Box aB1, aB2;
308   //
309   bSplit2 = Standard_False;
310   myRange1.Range(aT11, aT12);
311   myRange2.Range(aT21, aT22);
312   //
313   bIsClosed2 = IsClosed(myGeom2, aT21, aT22, myTol2, myRes2);
314   //
315   if (bIsClosed2) {
316     BndBuildBox(myCurve1, aT11, aT12, myTol1, aB1);
317     //
318     gp_Pnt aP = myGeom2->Value(aT21);
319     bIsClosed2 = !aB1.IsOut(aP);
320   }
321   //
322   if (!bIsClosed2) {
323     BndBuildBox(myCurve1, aT11, aT12, myTol1, aB1);
324     BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
325     FindSolutions(myRange1, aB1, myRange2, aB2, theRanges1, theRanges2);
326     return;
327   }
328   //
329   if (!CheckCoincidence(aT11, aT12, aT21, aT22, myTol, myRes1)) {
330     theRanges1.Append(myRange1);
331     theRanges2.Append(myRange2);
332     return;
333   }
334   //
335   Standard_Integer i, j, aNb1, aNb2;
336   IntTools_SequenceOfRanges aSegments1, aSegments2;
337   //
338   aNb1 = IsClosed(myGeom1, aT11, aT12, myTol1, myRes1) ? 2 : 1;
339   aNb2 = 2;
340   //
341   aNb1 = SplitRangeOnSegments(aT11, aT12, myRes1, aNb1, aSegments1);
342   aNb2 = SplitRangeOnSegments(aT21, aT22, myRes2, aNb2, aSegments2);
343   //
344   for (i = 1; i <= aNb1; ++i) {
345     const IntTools_Range& aR1 = aSegments1(i);
346     BndBuildBox(myCurve1, aR1.First(), aR1.Last(), myTol1, aB1);
347     for (j = 1; j <= aNb2; ++j) {
348       const IntTools_Range& aR2 = aSegments2(j);
349       BndBuildBox(myCurve2, aR2.First(), aR2.Last(), myTol2, aB2);
350       FindSolutions(aR1, aB1, aR2, aB2, theRanges1, theRanges2);
351     }
352   }
353   //
354   bSplit2 = aNb2 > 1;
355 }
356
357 //=======================================================================
358 //function : FindSolutions
359 //purpose  : 
360 //=======================================================================
361 void IntTools_EdgeEdge::FindSolutions(const IntTools_Range& theR1,
362                                       const Bnd_Box& theBox1,
363                                       const IntTools_Range& theR2,
364                                       const Bnd_Box& theBox2,
365                                       IntTools_SequenceOfRanges& theRanges1,
366                                       IntTools_SequenceOfRanges& theRanges2)
367 {
368   Standard_Boolean bOut, bStop, bThin;
369   Standard_Real aT11, aT12, aT21, aT22;
370   Standard_Real aTB11, aTB12, aTB21, aTB22;
371   Standard_Real aSmallStep1, aSmallStep2;
372   Standard_Integer iCom;
373   Bnd_Box aB1, aB2;
374   //
375   theR1.Range(aT11, aT12);
376   theR2.Range(aT21, aT22);
377   //
378   aB1 = theBox1;
379   aB2 = theBox2;
380   //
381   bThin = Standard_False;
382   bStop = Standard_False;
383   iCom  = 1;
384   //
385   do {
386     aTB11 = aT11;
387     aTB12 = aT12;
388     aTB21 = aT21;
389     aTB22 = aT22;
390     //
391     //1. Find parameters of the second edge in the box of first one
392     bOut = aB1.IsOut(aB2);
393     if (bOut) {
394       break;
395     }
396     //
397     bThin = ((aT12 - aT11) < myRes1) ||
398       (aB1.IsXThin(myTol) && aB1.IsYThin(myTol) && aB1.IsZThin(myTol));
399     //
400     bOut = !FindParameters(myCurve2, aTB21, aTB22, myTol2, myRes2, myPTol2, 
401                            myResCoeff2, aB1, aT21, aT22);
402     if (bOut || bThin) {
403       break;
404     }
405     //
406     //2. Build box for second edge and find parameters 
407     //   of the first one in that box
408     BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
409     bOut = aB1.IsOut(aB2);
410     if (bOut) {
411       break;
412     }
413     //
414     bThin = ((aT22 - aT21) < myRes2) ||
415       (aB2.IsXThin(myTol) && aB2.IsYThin(myTol) && aB2.IsZThin(myTol));
416     //
417     bOut = !FindParameters(myCurve1, aTB11, aTB12, myTol1, myRes1, myPTol1,
418                            myResCoeff1, aB2, aT11, aT12);
419     //
420     if (bOut || bThin) {
421       break;
422     }
423     //
424     //3. Check if it makes sense to continue
425     aSmallStep1 = (aTB12 - aTB11) / 250.;
426     aSmallStep2 = (aTB22 - aTB21) / 250.;
427     //
428     if (aSmallStep1 < myRes1) {
429       aSmallStep1 = myRes1;
430     }
431     if (aSmallStep2 < myRes2) {
432       aSmallStep2 = myRes2;
433     }
434     //
435     if (((aT11 - aTB11) < aSmallStep1) && ((aTB12 - aT12) < aSmallStep1) &&
436         ((aT21 - aTB21) < aSmallStep2) && ((aTB22 - aT22) < aSmallStep2)) {
437       bStop = Standard_True;
438     }
439     else
440       BndBuildBox (myCurve1, aT11, aT12, myTol1, aB1);
441
442   } while (!bStop);
443   //
444   if (bOut) {
445     //no intersection;
446     return;
447   }
448   //
449   if (!bThin) {
450     //check curves for coincidence on the ranges
451     iCom = CheckCoincidence(aT11, aT12, aT21, aT22, myTol, myRes1);
452     if (!iCom) {
453       bThin = Standard_True;
454     }
455   }
456   //
457   if (bThin) {
458     if (iCom != 0) {
459       //check intermediate points
460       Standard_Boolean bSol;
461       Standard_Real aT1;
462       gp_Pnt aP1;
463       GeomAPI_ProjectPointOnCurve aProjPC;
464       //
465       aT1 = (aT11 + aT12) * .5;
466       myGeom1->D0(aT1, aP1);
467       //
468       aProjPC.Init(myGeom2, aT21, aT22);
469       aProjPC.Perform(aP1);
470       //
471       if (aProjPC.NbPoints()) {
472         bSol = aProjPC.LowerDistance() <= myTol;
473       }
474       else {
475         Standard_Real aT2;
476         gp_Pnt aP2;
477         //
478         aT2 = (aT21 + aT22) * .5;
479         myGeom2->D0(aT2, aP2);
480         //
481         bSol = aP1.IsEqual(aP2, myTol);
482       }
483       //
484       if (!bSol) {
485         return;
486       }
487     }
488     //add common part
489     IntTools_Range aR1(aT11, aT12), aR2(aT21, aT22);
490     //
491     theRanges1.Append(aR1);
492     theRanges2.Append(aR2);
493     return;
494   }
495   //
496   if (!IsIntersection(aT11, aT12, aT21, aT22)) {
497     return;
498   }
499   //
500   //split ranges on segments and repeat
501   Standard_Integer i, aNb1;
502   IntTools_SequenceOfRanges aSegments1;
503   //
504   // Build box for first curve to compare
505   // the boxes of the splits with this one
506   BndBuildBox(myCurve1, aT11, aT12, myTol1, aB1);
507   const Standard_Real aB1SqExtent = aB1.SquareExtent();
508   //
509   IntTools_Range aR2(aT21, aT22);
510   BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
511   //
512   aNb1 = SplitRangeOnSegments(aT11, aT12, myRes1, 3, aSegments1);
513   for (i = 1; i <= aNb1; ++i) {
514     const IntTools_Range& aR1 = aSegments1(i);
515     BndBuildBox(myCurve1, aR1.First(), aR1.Last(), myTol1, aB1);
516     if (!aB1.IsOut(aB2) && (aNb1 == 1 || aB1.SquareExtent() < aB1SqExtent))
517       FindSolutions(aR1, aB1, aR2, aB2, theRanges1, theRanges2);
518   }
519 }
520
521 //=======================================================================
522 //function : FindParameters
523 //purpose  : 
524 //=======================================================================
525 Standard_Boolean IntTools_EdgeEdge::FindParameters(const BRepAdaptor_Curve& theBAC,
526                                                    const Standard_Real aT1, 
527                                                    const Standard_Real aT2,
528                                                    const Standard_Real theTol,
529                                                    const Standard_Real theRes,
530                                                    const Standard_Real thePTol,
531                                                    const Standard_Real theResCoeff,
532                                                    const Bnd_Box& theCBox,
533                                                    Standard_Real& aTB1, 
534                                                    Standard_Real& aTB2)
535 {
536   Standard_Boolean bRet;
537   Standard_Integer aC, i;
538   Standard_Real aCf, aDiff, aDt, aT, aTB, aTOut, aTIn;
539   Standard_Real aDist, aDistP;
540   gp_Pnt aP;
541   Bnd_Box aCBx;
542   //
543   bRet = Standard_False;
544   aCf = 0.6180339887498948482045868343656;// =0.5*(1.+sqrt(5.))/2.;
545   aCBx = theCBox;
546   aCBx.SetGap(aCBx.GetGap() + theTol);
547   //
548   const Handle(Geom_Curve)& aCurve = theBAC.Curve().Curve();
549   const GeomAbs_CurveType aCurveType = theBAC.GetType();
550   Standard_Real aMaxDt = (aT2 - aT1) * 0.01;
551   //
552   for (i = 0; i < 2; ++i) {
553     aTB = !i ? aT1 : aT2;
554     aT = !i ? aT2 : aTB1;
555     aC = !i ? 1 : -1;
556     aDt = theRes;
557     aDistP = 0.;
558     bRet = Standard_False;
559     Standard_Real k = 1;
560     //looking for the point on the edge which is in the box;
561     while (aC*(aT-aTB) >= 0) {
562       theBAC.D0(aTB, aP);
563       aDist = PointBoxDistance(theCBox, aP);
564       if (aDist > theTol) {
565         if (aDistP > 0.) {
566           Standard_Boolean toGrow = Standard_False;
567           if (Abs(aDistP - aDist) / aDistP < 0.1) {
568             aDt = Resolution(aCurve, aCurveType, theResCoeff, k*aDist);
569             if (aDt < aMaxDt)
570             {
571               toGrow = Standard_True;
572               k *= 2;
573             }
574           }
575           if (!toGrow) {
576             k = 1;
577             aDt = Resolution(aCurve, aCurveType, theResCoeff, aDist);
578           }
579         }
580         aTB += aC*aDt;
581       } else {
582         bRet = Standard_True;
583         break;
584       }
585       aDistP = aDist;
586     }
587     //
588     if (!bRet) {
589       if (!i) {
590         //edge is out of the box;
591         return bRet;
592       } else {
593         bRet = !bRet;
594         aTB = aTB1;
595         aDt = aT2 - aTB1;
596       }
597     }
598     //
599     aT = !i ? aT1 : aT2;
600     if (aTB != aT) {
601       //one point IN, one point OUT; looking for the bounding point;
602       aTIn = aTB;
603       aTOut = aTB - aC*aDt;
604       aDiff = aTIn - aTOut;
605       while (fabs(aDiff) > thePTol) {
606         aTB = aTOut + aDiff*aCf;
607         theBAC.D0(aTB, aP);
608         if (aCBx.IsOut(aP)) {
609           aTOut = aTB;
610         } else {
611           aTIn = aTB;
612         }
613         aDiff = aTIn - aTOut;
614       }
615     }
616     if (!i) {
617       aTB1 = aTB;
618     } else {
619       aTB2 = aTB;
620     }
621   }
622   return bRet;
623 }
624
625 //=======================================================================
626 //function : MergeSolutions
627 //purpose  : 
628 //=======================================================================
629 void IntTools_EdgeEdge::MergeSolutions(const IntTools_SequenceOfRanges& theRanges1, 
630                                        const IntTools_SequenceOfRanges& theRanges2,
631                                        const Standard_Boolean bSplit2)
632 {
633   Standard_Integer aNbCP = theRanges1.Length();
634   if (aNbCP == 0) {
635     return;
636   }
637   //
638   IntTools_Range aRi1, aRi2, aRj1, aRj2;
639   Standard_Boolean bCond;
640   Standard_Integer i, j;
641   TopAbs_ShapeEnum aType;
642   Standard_Real aT11, aT12, aT21, aT22;
643   Standard_Real aTi11, aTi12, aTi21, aTi22;
644   Standard_Real aTj11, aTj12, aTj21, aTj22;
645   Standard_Real aRes1, aRes2, dTR1, dTR2;
646   TColStd_MapOfInteger aMI;
647   //
648   aRes1 = Resolution(myCurve1.Curve().Curve(), 
649                      myCurve1.GetType(), myResCoeff1, myTol);
650   aRes2 = Resolution(myCurve2.Curve().Curve(), 
651                      myCurve2.GetType(), myResCoeff2, myTol);
652   //
653   myRange1.Range(aT11, aT12);
654   myRange2.Range(aT21, aT22);
655   dTR1 = 20*aRes1;
656   dTR2 = 20*aRes2;
657   aType = TopAbs_VERTEX;
658   //
659   for (i = 1; i <= aNbCP;) {
660     if (aMI.Contains(i)) {
661       ++i;
662       continue;
663     }
664     //
665     aRi1 = theRanges1(i);
666     aRi2 = theRanges2(i);
667     //
668     aRi1.Range(aTi11, aTi12);
669     aRi2.Range(aTi21, aTi22);
670     //
671     aMI.Add(i);
672     //
673     for (j = i+1; j <= aNbCP; ++j) {
674       if (aMI.Contains(j)) {
675         continue;
676       }
677       //
678       aRj1 = theRanges1(j);
679       aRj2 = theRanges2(j);
680       //
681       aRj1.Range(aTj11, aTj12);
682       aRj2.Range(aTj21, aTj22);
683       //
684       bCond = (fabs(aTi12 - aTj11) < dTR1) ||
685         (bSplit2 && (fabs(aTj12 - aTi11) < dTR1));
686       if (bCond && bSplit2) {
687         bCond = (fabs((Max(aTi22, aTj22) - Min(aTi21, aTj21)) - 
688                       ((aTi22 - aTi21) + (aTj22 - aTj21))) < dTR2);
689       }
690       //
691       if (bCond) {
692         aTi11 = Min(aTi11, aTj11);
693         aTi12 = Max(aTi12, aTj12);
694         aTi21 = Min(aTi21, aTj21);
695         aTi22 = Max(aTi22, aTj22);
696         aMI.Add(j);
697       }
698       else if (!bSplit2) {
699         i = j;
700         break;
701       }
702     }
703     //
704     if (((fabs(aT11 - aTi11) < myRes1) && (fabs(aT12 - aTi12) < myRes1)) ||
705         ((fabs(aT21 - aTi21) < myRes2) && (fabs(aT22 - aTi22) < myRes2))) {
706       aType = TopAbs_EDGE;
707       myCommonParts.Clear();
708     }
709     //
710     AddSolution(aTi11, aTi12, aTi21, aTi22, aType);
711     if (aType == TopAbs_EDGE) {
712       break;
713     }
714     //
715     if (bSplit2) {
716       ++i;
717     }
718   }
719 }
720
721 //=======================================================================
722 //function : AddSolution
723 //purpose  : 
724 //=======================================================================
725 void IntTools_EdgeEdge::AddSolution(const Standard_Real aT11,
726                                     const Standard_Real aT12,
727                                     const Standard_Real aT21,
728                                     const Standard_Real aT22,
729                                     const TopAbs_ShapeEnum theType)
730 {
731   IntTools_CommonPrt aCPart;
732   //
733   aCPart.SetType(theType);
734   if (!mySwap) {
735     aCPart.SetEdge1(myEdge1);
736     aCPart.SetEdge2(myEdge2);
737     aCPart.SetRange1(aT11, aT12);
738     aCPart.AppendRange2(aT21, aT22);
739   } else {
740     aCPart.SetEdge1(myEdge2);
741     aCPart.SetEdge2(myEdge1);
742     aCPart.SetRange1(aT21, aT22);
743     aCPart.AppendRange2(aT11, aT12);
744   }
745   //
746   if (theType == TopAbs_VERTEX) {
747     Standard_Real aT1, aT2;
748     //
749     FindBestSolution(aT11, aT12, aT21, aT22, aT1, aT2);
750     //
751     if (!mySwap) {
752       aCPart.SetVertexParameter1(aT1);
753       aCPart.SetVertexParameter2(aT2);
754     } else {
755       aCPart.SetVertexParameter1(aT2);
756       aCPart.SetVertexParameter2(aT1);
757     }
758   }
759   myCommonParts.Append(aCPart);
760 }
761
762 //=======================================================================
763 //function : FindBestSolution
764 //purpose  : 
765 //=======================================================================
766 void IntTools_EdgeEdge::FindBestSolution(const Standard_Real aT11,
767                                          const Standard_Real aT12,
768                                          const Standard_Real aT21,
769                                          const Standard_Real aT22,
770                                          Standard_Real& aT1,
771                                          Standard_Real& aT2)
772 {
773   Standard_Integer i, aNbS, iErr;
774   Standard_Real aDMin, aD, aRes1, aSolCriteria, aTouchCriteria;
775   Standard_Real aT1A, aT1B, aT1Min, aT2Min;
776   GeomAPI_ProjectPointOnCurve aProjPC;
777   IntTools_SequenceOfRanges aRanges;
778   //
779   aDMin = Precision::Infinite();
780   aSolCriteria   = 5.e-16;
781   aTouchCriteria = 5.e-13;
782   Standard_Boolean bTouch = Standard_False;
783   Standard_Boolean bTouchConfirm = Standard_False;
784   //
785   aRes1 = Resolution(myCurve1.Curve().Curve(), 
786                      myCurve1.GetType(), myResCoeff1, myTol);
787   aNbS = 10;
788   aNbS = SplitRangeOnSegments(aT11, aT12, 3*aRes1, aNbS, aRanges);
789   //
790   aProjPC.Init(myGeom2, aT21, aT22);
791   //
792   Standard_Real aT11Touch = aT11, aT12Touch = aT12;
793   Standard_Real aT21Touch = aT21, aT22Touch = aT22;
794   Standard_Boolean isSolFound = Standard_False;
795   for (i = 1; i <= aNbS; ++i) {
796     const IntTools_Range& aR1 = aRanges(i);
797     aR1.Range(aT1A, aT1B);
798     //
799     aD = myTol;
800     iErr = FindDistPC(aT1A, aT1B, myGeom1, aSolCriteria, myPTol1,
801                       aProjPC, aD, aT1Min, aT2Min, Standard_False);
802     if (iErr != 1) {
803       if (aD < aDMin) {
804         aT1 = aT1Min;
805         aT2 = aT2Min;
806         aDMin = aD;
807         isSolFound = Standard_True;
808       }
809       //
810       if (aD < aTouchCriteria) {
811         if (bTouch) {
812           aT12Touch = aT1Min;
813           aT22Touch = aT2Min;
814           bTouchConfirm = Standard_True;
815         }
816         else {
817           aT11Touch = aT1Min;
818           aT21Touch = aT2Min;
819           bTouch = Standard_True;
820         }
821       }
822     }
823   }
824   if (!isSolFound || bTouchConfirm)
825   {
826     aT1 = (aT11Touch + aT12Touch) * 0.5;
827     iErr = DistPC(aT1, myGeom1, aSolCriteria, aProjPC, aD, aT2, -1);
828     if (iErr == 1) {
829       aT2 = (aT21Touch + aT22Touch) * 0.5;
830     }
831   }
832 }
833
834 //=======================================================================
835 //function : ComputeLineLine
836 //purpose  : 
837 //=======================================================================
838 void IntTools_EdgeEdge::ComputeLineLine()
839 {
840   Standard_Real aTol = myTol * myTol;
841
842   gp_Lin aL1 = myCurve1.Line();
843   gp_Lin aL2 = myCurve2.Line();
844
845   gp_Dir aD1 = aL1.Direction();
846   gp_Dir aD2 = aL2.Direction();
847
848   Standard_Real anAngle = aD1.Angle (aD2);
849   Standard_Boolean IsCoincide = anAngle < Precision::Angular();
850   if (IsCoincide)
851   {
852     if (aL1.SquareDistance (aL2.Location()) > aTol)
853       return;
854   }
855
856   Standard_Real aT11, aT12, aT21, aT22;
857   myRange1.Range (aT11, aT12);
858   myRange2.Range (aT21, aT22);
859
860   gp_Pnt aP11 = ElCLib::Value (aT11, aL1);
861   gp_Pnt aP12 = ElCLib::Value (aT12, aL1);
862
863   if (!IsCoincide)
864   {
865     gp_Pnt O2 (aL2.Location());
866     if (!Precision::IsInfinite (aT21) && !Precision::IsInfinite (aT22))
867       O2 = ElCLib::Value ((aT21 + aT22) / 2., aL2);
868
869     gp_Vec aVec1 = gp_Vec (O2, aP11).Crossed (aD2);
870     gp_Vec aVec2 = gp_Vec (O2, aP12).Crossed (aD2);
871
872     Standard_Real aSqDist1 = aVec1.SquareMagnitude();
873     Standard_Real aSqDist2 = aVec2.SquareMagnitude();
874
875     IsCoincide = (aSqDist1 <= aTol && aSqDist2 <= aTol);
876
877     if (!IsCoincide && aVec1.Dot (aVec2) > 0)
878       // the lines do not intersect
879       return;
880   }
881
882   IntTools_CommonPrt aCommonPrt;
883   aCommonPrt.SetEdge1 (myEdge1);
884   aCommonPrt.SetEdge2 (myEdge2);
885
886   if (IsCoincide)
887   {
888     Standard_Real t21 = ElCLib::Parameter (aL2, aP11);
889     Standard_Real t22 = ElCLib::Parameter (aL2, aP12);
890
891     if ((t21 > aT22 && t22 > aT22) || (t21 < aT21 && t22 < aT21))
892       // projections are out of range
893       return;
894
895     if (t21 > t22)
896       std::swap (t21, t22);
897
898     if (t21 >= aT21)
899     {
900       if (t22 <= aT22)
901       {
902         aCommonPrt.SetRange1 (aT11, aT12);
903         aCommonPrt.SetAllNullFlag (Standard_True);
904         aCommonPrt.AppendRange2 (t21, t22);
905       }
906       else
907       {
908         aCommonPrt.SetRange1 (aT11, aT12 - (t22 - aT22));
909         aCommonPrt.AppendRange2 (t21, aT22);
910       }
911     }
912     else
913     {
914       aCommonPrt.SetRange1 (aT11 + (aT21 - t21), aT12);
915       aCommonPrt.AppendRange2 (aT21, t22);
916     }
917     aCommonPrt.SetType (TopAbs_EDGE);
918     myCommonParts.Append (aCommonPrt);
919     return;
920   }
921
922
923   gp_Vec O1O2 (aL1.Location(), aL2.Location());
924   gp_XYZ aCross = aD1.XYZ().Crossed (aD2.XYZ());
925   Standard_Real aDistLL = O1O2.Dot (gp_Vec (aCross.Normalized()));
926   if (Abs (aDistLL) > myTol)
927     return;
928
929   {
930     // Fast check that no intersection needs to be added
931     for (TopoDS_Iterator it1 (myEdge1); it1.More(); it1.Next())
932     {
933       for (TopoDS_Iterator it2 (myEdge2); it2.More(); it2.Next())
934       {
935         if (it1.Value().IsSame (it2.Value()))
936           return;
937       }
938     }
939   }
940
941   Standard_Real aSqSin = aCross.SquareModulus();
942   Standard_Real aT2 = (aD1.XYZ() * (O1O2.Dot (aD1)) - (O1O2.XYZ())).Dot (aD2.XYZ());
943   aT2 /= aSqSin;
944
945   if (aT2 < aT21 || aT2 > aT22)
946     // out of range
947     return;
948
949   gp_Pnt aP2 = ElCLib::Value (aT2, aL2);
950   Standard_Real aT1 = gp_Vec (aL1.Location(), aP2).Dot (aD1);
951
952   if (aT1 < aT11 || aT1 > aT12)
953     // out of range
954     return;
955
956   gp_Pnt aP1 = ElCLib::Value (aT1, aL1);
957   Standard_Real aDist = aP1.SquareDistance (aP2);
958
959   if (aDist > aTol)
960     // no intersection
961     return;
962
963   // compute correct range on the edges
964   Standard_Real aDt1 = IntTools_Tools::ComputeIntRange (myTol1, myTol2, anAngle);
965   Standard_Real aDt2 = IntTools_Tools::ComputeIntRange (myTol2, myTol1, anAngle);
966
967   aCommonPrt.SetRange1 (aT1 - aDt1, aT1 + aDt1);
968   aCommonPrt.AppendRange2 (aT2 - aDt2, aT2 + aDt2);
969   aCommonPrt.SetType (TopAbs_VERTEX);
970   aCommonPrt.SetVertexParameter1 (aT1);
971   aCommonPrt.SetVertexParameter2 (aT2);
972   myCommonParts.Append (aCommonPrt);
973 }
974
975 //=======================================================================
976 //function : IsIntersection
977 //purpose  : 
978 //=======================================================================
979 Standard_Boolean IntTools_EdgeEdge::IsIntersection(const Standard_Real aT11,
980                                                    const Standard_Real aT12,
981                                                    const Standard_Real aT21,
982                                                    const Standard_Real aT22)
983 {
984   Standard_Boolean bRet;
985   gp_Pnt aP11, aP12, aP21, aP22;
986   gp_Vec aV11, aV12, aV21, aV22;
987   Standard_Real aD11_21, aD11_22, aD12_21, aD12_22, aCriteria, aCoef;
988   Standard_Boolean bSmall_11_21, bSmall_11_22, bSmall_12_21, bSmall_12_22;
989   //
990   bRet = Standard_True;
991   aCoef = 1.e+5;
992   if (((aT12 - aT11) > aCoef*myRes1) && ((aT22 - aT21) > aCoef*myRes2)) {
993     aCoef = 5000;
994   } else {
995     Standard_Real aTRMin = Min((aT12 - aT11)/myRes1, (aT22 - aT21)/myRes2);
996     aCoef = aTRMin / 100.;
997     if (aCoef < 1.) {
998       aCoef = 1.;
999     }
1000   }
1001   aCriteria = aCoef * myTol;
1002   aCriteria *= aCriteria;
1003   //
1004   myGeom1->D1(aT11, aP11, aV11);
1005   myGeom1->D1(aT12, aP12, aV12);
1006   myGeom2->D1(aT21, aP21, aV21);
1007   myGeom2->D1(aT22, aP22, aV22);
1008   //
1009   aD11_21 = aP11.SquareDistance(aP21);
1010   aD11_22 = aP11.SquareDistance(aP22);
1011   aD12_21 = aP12.SquareDistance(aP21);
1012   aD12_22 = aP12.SquareDistance(aP22);
1013   //
1014   bSmall_11_21 = aD11_21 < aCriteria;
1015   bSmall_11_22 = aD11_22 < aCriteria;
1016   bSmall_12_21 = aD12_21 < aCriteria;
1017   bSmall_12_22 = aD12_22 < aCriteria;
1018   //
1019   if ((bSmall_11_21 && bSmall_12_22) ||
1020       (bSmall_11_22 && bSmall_12_21)) {
1021     if (aCoef == 1.) {
1022       return bRet;
1023     }
1024     //
1025     Standard_Real anAngleCriteria;
1026     Standard_Real anAngle1 = 0.0,
1027                   anAngle2 = 0.0;
1028     //
1029     anAngleCriteria = 5.e-3;
1030     if (aV11.SquareMagnitude() > Precision::SquareConfusion() &&
1031         aV12.SquareMagnitude() > Precision::SquareConfusion() &&
1032         aV21.SquareMagnitude() > Precision::SquareConfusion() &&
1033         aV22.SquareMagnitude() > Precision::SquareConfusion() )
1034     {
1035       if (bSmall_11_21 && bSmall_12_22) {
1036         anAngle1 = aV11.Angle(aV21);
1037         anAngle2 = aV12.Angle(aV22);
1038       } else {
1039         anAngle1 = aV11.Angle(aV22);
1040         anAngle2 = aV12.Angle(aV21);
1041       }
1042     }
1043     //
1044     if (((anAngle1 < anAngleCriteria) || ((M_PI - anAngle1) < anAngleCriteria)) ||
1045         ((anAngle2 < anAngleCriteria) || ((M_PI - anAngle2) < anAngleCriteria))) {
1046       GeomAPI_ProjectPointOnCurve aProjPC;
1047       Standard_Integer iErr;
1048       Standard_Real aD, aT1Min, aT2Min;
1049       //
1050       aD = Precision::Infinite();
1051       aProjPC.Init(myGeom2, aT21, aT22);
1052       iErr = FindDistPC(aT11, aT12, myGeom1, myTol, myRes1, 
1053                         aProjPC, aD, aT1Min, aT2Min, Standard_False);
1054       bRet = (iErr == 2);
1055     }
1056   }
1057   return bRet;
1058 }
1059
1060 //=======================================================================
1061 //function : CheckCoincidence
1062 //purpose  : 
1063 //=======================================================================
1064 Standard_Integer IntTools_EdgeEdge::CheckCoincidence(const Standard_Real aT11,
1065                                                      const Standard_Real aT12,
1066                                                      const Standard_Real aT21,
1067                                                      const Standard_Real aT22,
1068                                                      const Standard_Real theCriteria,
1069                                                      const Standard_Real theCurveRes1)
1070 {
1071   Standard_Integer iErr, aNb, aNb1, i;
1072   Standard_Real aT1A, aT1B, aT1max, aT2max, aDmax;
1073   GeomAPI_ProjectPointOnCurve aProjPC;
1074   IntTools_SequenceOfRanges aRanges;
1075   //
1076   iErr  = 0;
1077   aDmax = -1.;
1078   aProjPC.Init(myGeom2, aT21, aT22);
1079   //
1080   // 1. Express evaluation
1081   aNb = 10; // Number of intervals on the curve #1
1082   aNb1 = SplitRangeOnSegments(aT11, aT12, theCurveRes1, aNb, aRanges);
1083   for (i = 1; i < aNb1; ++i) {
1084     const IntTools_Range& aR1 = aRanges(i);
1085     aR1.Range(aT1A, aT1B);
1086     //
1087     iErr = DistPC(aT1B, myGeom1, theCriteria, aProjPC, aDmax, aT2max);
1088     if (iErr) {
1089       return iErr;
1090     }
1091   }
1092   //
1093   // if the ranges in aRanges are less than theCurveRes1,
1094   // there is no need to do step 2 (deep evaluation)
1095   if (aNb1 < aNb) {
1096     return iErr;
1097   }
1098   //
1099   // 2. Deep evaluation
1100   for (i = 2; i < aNb1; ++i) {
1101     const IntTools_Range& aR1 = aRanges(i);
1102     aR1.Range(aT1A, aT1B);
1103     //
1104     iErr = FindDistPC(aT1A, aT1B, myGeom1, theCriteria, theCurveRes1, 
1105                       aProjPC, aDmax, aT1max, aT2max);
1106     if (iErr) {
1107       return iErr;
1108     }
1109   }
1110   // Possible values:
1111   // iErr == 0 - the patches are coincided
1112   // iErr == 1 - a point from aC1 can not be projected on aC2
1113   // iErr == 2 - the distance is too big
1114   return iErr;
1115 }
1116
1117 //=======================================================================
1118 //function : FindDistPC
1119 //purpose  : 
1120 //=======================================================================
1121 Standard_Integer FindDistPC(const Standard_Real aT1A, 
1122                             const Standard_Real aT1B,
1123                             const Handle(Geom_Curve)& theC1,
1124                             const Standard_Real theCriteria,
1125                             const Standard_Real theEps,
1126                             GeomAPI_ProjectPointOnCurve& theProjPC,
1127                             Standard_Real& aDmax, 
1128                             Standard_Real& aT1max,
1129                             Standard_Real& aT2max,
1130                             const Standard_Boolean bMaxDist) 
1131 {
1132   Standard_Integer iErr, iC;
1133   Standard_Real aGS, aXP, aA, aB, aXL, aYP, aYL, aT2P, aT2L;
1134   //
1135   iC = bMaxDist ? 1 : -1;
1136   iErr = 0;
1137   aT1max = aT2max = 0.; // silence GCC warning
1138   //
1139   aGS = 0.6180339887498948482045868343656;// =0.5*(1.+sqrt(5.))-1.;
1140   aA = aT1A;
1141   aB = aT1B;
1142   //
1143   // check bounds
1144   iErr = DistPC(aA, theC1, theCriteria, theProjPC, 
1145                 aYP, aT2P, aDmax, aT1max, aT2max, iC);
1146   if (iErr == 2) {
1147     return iErr;
1148   }
1149   //
1150   iErr = DistPC(aB, theC1, theCriteria, theProjPC, 
1151                 aYL, aT2L, aDmax, aT1max, aT2max, iC);
1152   if (iErr == 2) {
1153     return iErr;
1154   }
1155   //
1156   aXP = aA + (aB - aA)*aGS;
1157   aXL = aB - (aB - aA)*aGS;
1158   //
1159   iErr = DistPC(aXP, theC1, theCriteria, theProjPC, 
1160                 aYP, aT2P, aDmax, aT1max, aT2max, iC);
1161   if (iErr) {
1162     return iErr;
1163   }
1164   //
1165   iErr = DistPC(aXL, theC1, theCriteria, theProjPC, 
1166                 aYL, aT2L, aDmax, aT1max, aT2max, iC);
1167   if (iErr) {
1168     return iErr;
1169   }
1170   //
1171   Standard_Real anEps = Max(theEps, Epsilon(Max(Abs(aA), Abs(aB))) * 10.);
1172   for (;;) {
1173     if (iC*(aYP - aYL) > 0) {
1174       aA = aXL;
1175       aXL = aXP;
1176       aYL = aYP;
1177       aXP = aA + (aB - aA)*aGS;
1178       iErr = DistPC(aXP, theC1, theCriteria, theProjPC, 
1179                     aYP, aT2P, aDmax, aT1max, aT2max, iC);
1180     }
1181     else {
1182       aB = aXP;
1183       aXP = aXL;
1184       aYP = aYL;
1185       aXL = aB - (aB - aA)*aGS;
1186       iErr = DistPC(aXL, theC1, theCriteria, theProjPC, 
1187                     aYL, aT2L, aDmax, aT1max, aT2max, iC);
1188     }
1189     //
1190     if (iErr) {
1191       if ((iErr == 2) && !bMaxDist) {
1192         aXP = (aA + aB) * 0.5;
1193         DistPC(aXP, theC1, theCriteria, theProjPC, 
1194                aYP, aT2P, aDmax, aT1max, aT2max, iC);
1195       }
1196       return iErr;
1197     }
1198     //
1199     if ((aB - aA) < anEps) {
1200       break;
1201     }
1202   }// for (;;) {
1203   //
1204   return iErr;
1205 }
1206 //=======================================================================
1207 //function : DistPC
1208 //purpose  : 
1209 //=======================================================================
1210 Standard_Integer DistPC(const Standard_Real aT1, 
1211                         const Handle(Geom_Curve)& theC1,
1212                         const Standard_Real theCriteria,
1213                         GeomAPI_ProjectPointOnCurve& theProjPC, 
1214                         Standard_Real& aD, 
1215                         Standard_Real& aT2,
1216                         Standard_Real& aDmax,
1217                         Standard_Real& aT1max,
1218                         Standard_Real& aT2max,
1219                         const Standard_Integer iC)
1220 {
1221   Standard_Integer iErr;
1222   //
1223   iErr = DistPC(aT1, theC1, theCriteria, theProjPC, aD, aT2, iC);
1224   if (iErr == 1) {
1225     return iErr;
1226   }
1227   //
1228   if (iC*(aD - aDmax) > 0) {
1229     aDmax = aD;
1230     aT1max = aT1;
1231     aT2max = aT2;
1232   }
1233   //
1234   return iErr;
1235 }
1236 //=======================================================================
1237 //function : DistPC
1238 //purpose  : 
1239 //=======================================================================
1240 Standard_Integer DistPC(const Standard_Real aT1, 
1241                         const Handle(Geom_Curve)& theC1,
1242                         const Standard_Real theCriteria, 
1243                         GeomAPI_ProjectPointOnCurve& theProjPC,
1244                         Standard_Real& aD, 
1245                         Standard_Real& aT2,
1246                         const Standard_Integer iC) 
1247 {
1248   Standard_Integer iErr, aNbP2;
1249   gp_Pnt aP1;
1250   //
1251   iErr = 0;
1252   theC1->D0(aT1, aP1);
1253   //
1254   theProjPC.Perform(aP1);
1255   aNbP2 = theProjPC.NbPoints();
1256   if (!aNbP2) {
1257     iErr = 1;// the point from aC1 can not be projected on aC2
1258     return iErr;
1259   }
1260   //
1261   aD  = theProjPC.LowerDistance();
1262   aT2 = theProjPC.LowerDistanceParameter();
1263   if (iC*(aD - theCriteria) > 0) {
1264     iErr = 2;// the distance is too big or small
1265   }
1266   //
1267   return iErr;
1268 }
1269
1270 //=======================================================================
1271 //function : SplitRangeOnSegments
1272 //purpose  : 
1273 //=======================================================================
1274 Standard_Integer SplitRangeOnSegments(const Standard_Real aT1, 
1275                                       const Standard_Real aT2,
1276                                       const Standard_Real theResolution,
1277                                       const Standard_Integer theNbSeg,
1278                                       IntTools_SequenceOfRanges& theSegments)
1279 {
1280   Standard_Real aDiff = aT2 - aT1;
1281   if (aDiff < theResolution || theNbSeg == 1) {
1282     theSegments.Append(IntTools_Range(aT1, aT2));
1283     return 1;
1284   }
1285   //
1286   Standard_Real aDt, aT1x, aT2x, aSeg;
1287   Standard_Integer aNbSegments, i;
1288   //
1289   aNbSegments = theNbSeg;
1290   aDt = aDiff / aNbSegments;
1291   if (aDt < theResolution) {
1292     aSeg = aDiff / theResolution;
1293     aNbSegments = Standard_Integer(aSeg) + 1;
1294     aDt = aDiff / aNbSegments;
1295   }
1296   //
1297   aT1x = aT1;
1298   for (i = 1; i < aNbSegments; ++i) {
1299     aT2x = aT1x + aDt;
1300     //
1301     IntTools_Range aR(aT1x, aT2x);
1302     theSegments.Append(aR);
1303     //
1304     aT1x = aT2x;
1305   }
1306   //
1307   IntTools_Range aR(aT1x, aT2);
1308   theSegments.Append(aR);
1309   //
1310   return aNbSegments;
1311 }
1312
1313 //=======================================================================
1314 //function : BndBuildBox
1315 //purpose  : 
1316 //=======================================================================
1317 void BndBuildBox(const BRepAdaptor_Curve& theBAC,
1318                  const Standard_Real aT1,
1319                  const Standard_Real aT2,
1320                  const Standard_Real theTol,
1321                  Bnd_Box& theBox)
1322 {
1323   Bnd_Box aB;
1324   BndLib_Add3dCurve::Add(theBAC, aT1, aT2, theTol, aB);
1325   theBox = aB;
1326 }
1327
1328 //=======================================================================
1329 //function : PointBoxDistance
1330 //purpose  : 
1331 //=======================================================================
1332 Standard_Real PointBoxDistance(const Bnd_Box& aB,
1333                                const gp_Pnt& aP)
1334 {
1335   Standard_Real aPCoord[3];
1336   Standard_Real aBMinCoord[3], aBMaxCoord[3];
1337   Standard_Real aDist, aR1, aR2;
1338   Standard_Integer i;
1339   //
1340   aP.Coord(aPCoord[0], aPCoord[1], aPCoord[2]);
1341   aB.Get(aBMinCoord[0], aBMinCoord[1], aBMinCoord[2], 
1342          aBMaxCoord[0], aBMaxCoord[1], aBMaxCoord[2]);
1343   //
1344   aDist = 0.;
1345   for (i = 0; i < 3; ++i) {
1346     aR1 = aBMinCoord[i] - aPCoord[i];
1347     if (aR1 > 0.) {
1348       aDist += aR1*aR1;
1349       continue;
1350     }
1351     //
1352     aR2 = aPCoord[i] - aBMaxCoord[i];
1353     if (aR2 > 0.) {
1354       aDist += aR2*aR2;
1355     }
1356   }
1357   //
1358   aDist = Sqrt(aDist);
1359   return aDist;
1360 }
1361
1362 //=======================================================================
1363 //function : TypeToInteger
1364 //purpose  : 
1365 //=======================================================================
1366 Standard_Integer TypeToInteger(const GeomAbs_CurveType theCType)
1367 {
1368   Standard_Integer iRet;
1369   //
1370   switch(theCType) {
1371   case GeomAbs_Line:
1372     iRet=0;
1373     break;
1374   case GeomAbs_Hyperbola:
1375   case GeomAbs_Parabola:
1376     iRet=1;
1377     break;
1378   case GeomAbs_Circle:
1379   case GeomAbs_Ellipse:
1380     iRet=2;
1381     break;
1382   case GeomAbs_BezierCurve:
1383   case GeomAbs_BSplineCurve:
1384     iRet=3;
1385     break;
1386   default:
1387     iRet=4;
1388     break;
1389   }
1390   return iRet;
1391 }
1392
1393 //=======================================================================
1394 //function : ResolutionCoeff
1395 //purpose  : 
1396 //=======================================================================
1397 Standard_Real ResolutionCoeff(const BRepAdaptor_Curve& theBAC,
1398                               const IntTools_Range& theRange)
1399 {
1400   Standard_Real aResCoeff = 0.;
1401   //
1402   const Handle(Geom_Curve)& aCurve = theBAC.Curve().Curve();
1403   const GeomAbs_CurveType aCurveType = theBAC.GetType();
1404   //
1405   switch (aCurveType) {
1406   case GeomAbs_Circle :
1407     aResCoeff = 1. / (2 * Handle(Geom_Circle)::DownCast (aCurve)->Circ().Radius());
1408     break;
1409   case GeomAbs_Ellipse :
1410     aResCoeff =  1. / Handle(Geom_Ellipse)::DownCast (aCurve)->MajorRadius();
1411     break;
1412   case GeomAbs_OffsetCurve : {
1413     const Handle(Geom_OffsetCurve)& anOffsetCurve = Handle(Geom_OffsetCurve)::DownCast(aCurve);
1414     const Handle(Geom_Curve)& aBasisCurve = anOffsetCurve->BasisCurve();
1415     GeomAdaptor_Curve aGBasisCurve(aBasisCurve);
1416     const GeomAbs_CurveType aBCType = aGBasisCurve.GetType();
1417     if (aBCType == GeomAbs_Line) {
1418       break;
1419     }
1420     else if (aBCType == GeomAbs_Circle) {
1421       aResCoeff = 1. / (2 * (anOffsetCurve->Offset() + aGBasisCurve.Circle().Radius()));
1422       break;
1423     }
1424     else if (aBCType == GeomAbs_Ellipse) {
1425       aResCoeff = 1. / (anOffsetCurve->Offset() + aGBasisCurve.Ellipse().MajorRadius());
1426       break;
1427     }
1428   }
1429   Standard_FALLTHROUGH
1430   case GeomAbs_Hyperbola :
1431   case GeomAbs_Parabola : 
1432   case GeomAbs_OtherCurve :{
1433     Standard_Real k, kMin, aDist, aDt, aT1, aT2, aT;
1434     Standard_Integer aNbP, i;
1435     gp_Pnt aP1, aP2;
1436     //
1437     aNbP = 30;
1438     theRange.Range(aT1, aT2);
1439     aDt = (aT2 - aT1) / aNbP;
1440     aT = aT1;
1441     kMin = 10.;
1442     //
1443     theBAC.D0(aT1, aP1);
1444     for (i = 1; i <= aNbP; ++i) {
1445       aT += aDt;
1446       theBAC.D0(aT, aP2);
1447       aDist = aP1.Distance(aP2);
1448       k = aDt / aDist;
1449       if (k < kMin) {
1450         kMin = k;
1451       }
1452       aP1 = aP2;
1453     }
1454     //
1455     aResCoeff = kMin;
1456     break;
1457   }
1458   default:
1459     break;
1460   }
1461   //
1462   return aResCoeff;
1463 }
1464
1465 //=======================================================================
1466 //function : Resolution
1467 //purpose  : 
1468 //=======================================================================
1469 Standard_Real Resolution(const Handle(Geom_Curve)& theCurve,
1470                          const GeomAbs_CurveType theCurveType,
1471                          const Standard_Real theResCoeff,
1472                          const Standard_Real theR3D)
1473 {
1474   Standard_Real aRes;
1475   //
1476   switch (theCurveType) {
1477   case GeomAbs_Line :
1478     aRes = theR3D;
1479     break;
1480   case GeomAbs_Circle: {
1481     Standard_Real aDt = theResCoeff * theR3D;
1482     aRes = (aDt <= 1.) ? 2*ASin(aDt) : 2*M_PI;
1483     break;
1484   }
1485   case GeomAbs_BezierCurve:
1486     Handle(Geom_BezierCurve)::DownCast (theCurve)->Resolution(theR3D, aRes);
1487     break;
1488   case GeomAbs_BSplineCurve:
1489     Handle(Geom_BSplineCurve)::DownCast (theCurve)->Resolution(theR3D, aRes);
1490     break;
1491   case GeomAbs_OffsetCurve: {
1492     const Handle(Geom_Curve)& aBasisCurve = 
1493       Handle(Geom_OffsetCurve)::DownCast(theCurve)->BasisCurve();
1494     const GeomAbs_CurveType aBCType = GeomAdaptor_Curve(aBasisCurve).GetType();
1495     if (aBCType == GeomAbs_Line) {
1496       aRes = theR3D;
1497       break;
1498     }
1499     else if (aBCType == GeomAbs_Circle) {
1500       Standard_Real aDt = theResCoeff * theR3D;
1501       aRes = (aDt <= 1.) ? 2*ASin(aDt) : 2*M_PI;
1502       break;
1503     }
1504   }
1505   Standard_FALLTHROUGH
1506   default:
1507     aRes = theResCoeff * theR3D;
1508     break;
1509   }
1510   //
1511   return aRes;
1512 }
1513
1514 //=======================================================================
1515 //function : CurveDeflection
1516 //purpose  : 
1517 //=======================================================================
1518 Standard_Real CurveDeflection(const BRepAdaptor_Curve& theBAC,
1519                               const IntTools_Range& theRange)
1520 {
1521   Standard_Real aDt, aT, aT1, aT2, aDefl;
1522   Standard_Integer i, aNbP;
1523   gp_Vec aV1, aV2;
1524   gp_Pnt aP;
1525   //
1526   aDefl = 0;
1527   aNbP = 10;
1528   theRange.Range(aT1, aT2);
1529   aDt = (aT2 - aT1) / aNbP;
1530   aT = aT1;
1531   //
1532   theBAC.D1(aT1, aP, aV1);
1533   for (i = 1; i <= aNbP; ++i) {
1534     aT += aDt;
1535     theBAC.D1(aT, aP, aV2);
1536     if (aV1.Magnitude() > gp::Resolution() &&
1537         aV2.Magnitude() > gp::Resolution()) {
1538       gp_Dir aD1(aV1), aD2(aV2);
1539       aDefl += aD1.Angle(aD2);
1540     }
1541     aV1 = aV2;
1542   }
1543   //
1544   return aDefl;
1545 }
1546
1547 //=======================================================================
1548 //function : IsClosed
1549 //purpose  : 
1550 //=======================================================================
1551 Standard_Boolean IsClosed(const Handle(Geom_Curve)& theCurve,
1552                           const Standard_Real aT1,
1553                           const Standard_Real aT2,
1554                           const Standard_Real theTol,
1555                           const Standard_Real theRes)
1556 {
1557   if (Abs(aT1 - aT2) < theRes)
1558   {
1559     return Standard_False;
1560   }
1561
1562   gp_Pnt aP1, aP2;
1563   theCurve->D0(aT1, aP1);
1564   theCurve->D0(aT2, aP2);
1565   //
1566   Standard_Real aD = aP1.Distance(aP2);
1567   return aD < theTol;
1568 }