44bc9a0884c8cd3ea81881247be8a35e0bb0e630
[occt.git] / src / IntPatch / IntPatch_ALineToWLine.cxx
1 // Created on: 1993-11-26
2 // Created by: Modelistation
3 // Copyright (c) 1993-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #include <IntPatch_ALineToWLine.hxx>
18
19 #include <Adaptor3d_HSurface.hxx>
20 #include <ElSLib.hxx>
21 #include <IntPatch_ALine.hxx>
22 #include <IntPatch_Point.hxx>
23 #include <IntPatch_SpecialPoints.hxx>
24 #include <IntPatch_WLine.hxx>
25 #include <IntSurf.hxx>
26 #include <IntSurf_LineOn2S.hxx>
27
28 //=======================================================================
29 //function : AddPointIntoLine
30 //purpose  : 
31 //=======================================================================
32 static inline void AddPointIntoLine(Handle(IntSurf_LineOn2S) theLine,
33                                     const Standard_Real* const theArrPeriods,
34                                     IntSurf_PntOn2S &thePoint,
35                                     IntPatch_Point* theVertex = 0)
36 {
37    if(theLine->NbPoints() > 0)
38    {
39      if(thePoint.IsSame(theLine->Value(theLine->NbPoints()), Precision::Confusion()))
40        return;
41
42      IntPatch_SpecialPoints::AdjustPointAndVertex(theLine->Value(theLine->NbPoints()),
43                                                     theArrPeriods, thePoint, theVertex);
44    }
45
46    theLine->Add(thePoint);
47 }
48
49 //=======================================================================
50 //function : AddVertexPoint
51 //purpose  : Extracts IntSurf_PntOn2S from theVertex and adds result in theLine.
52 //=======================================================================
53 static void AddVertexPoint(Handle(IntSurf_LineOn2S)& theLine,
54                                        IntPatch_Point &theVertex,
55                                        const Standard_Real* const theArrPeriods)
56 {
57   IntSurf_PntOn2S anApexPoint = theVertex.PntOn2S();
58   AddPointIntoLine(theLine, theArrPeriods, anApexPoint, &theVertex);
59 }
60
61 //=======================================================================
62 //function : IsPoleOrSeam
63 //purpose  : Processes theVertex depending on its type
64 //            (pole/apex/point on boundary etc.) and adds it in theLine.
65 //           thePIsoRef is the reference point using in case when the
66 //            value of correspond parameter cannot be precise.
67 //           theSingularSurfaceID contains the ID of surface with
68 //            special point (0 - none, 1 - theS1, 2 - theS2)
69 //=======================================================================
70 static IntPatch_SpecPntType IsPoleOrSeam(const Handle(Adaptor3d_HSurface)& theS1,
71                                          const Handle(Adaptor3d_HSurface)& theS2,
72                                          const IntSurf_PntOn2S& thePIsoRef,
73                                          Handle(IntSurf_LineOn2S)& theLine,
74                                          IntPatch_Point &theVertex,
75                                          const Standard_Real theArrPeriods[4],
76                                          const Standard_Real theTol3d,
77                                          Standard_Integer& theSingularSurfaceID)
78 {
79   theSingularSurfaceID = 0;
80
81   for(Standard_Integer i = 0; i < 2; i++)
82   {
83     const Standard_Boolean isReversed = (i > 0);
84     const GeomAbs_SurfaceType aType = isReversed? theS2->GetType() : theS1->GetType();
85
86     IntPatch_SpecPntType anAddedPType = IntPatch_SPntNone;
87     IntSurf_PntOn2S anApexPoint;
88
89     switch(aType)
90     {
91     case GeomAbs_Sphere:
92     case GeomAbs_Cone:
93       {
94         if(IntPatch_SpecialPoints::
95               AddSingularPole((isReversed? theS2 : theS1), (isReversed? theS1 : theS2),
96                               thePIsoRef, theVertex, anApexPoint,
97                               isReversed, Standard_True))
98         {
99           anAddedPType = IntPatch_SPntPole;
100           break;
101         }
102       }
103       Standard_FALLTHROUGH
104     case GeomAbs_Torus:
105       if(aType == GeomAbs_Torus)
106       {
107         if(IntPatch_SpecialPoints::
108             AddCrossUVIsoPoint((isReversed? theS2 : theS1), (isReversed? theS1 : theS2),
109                                   thePIsoRef, theTol3d,
110                                   anApexPoint, isReversed))
111         {
112           anAddedPType = IntPatch_SPntSeamUV;
113           break;
114         }
115       }
116       Standard_FALLTHROUGH
117     case GeomAbs_Cylinder:
118       theSingularSurfaceID = i + 1;
119       AddVertexPoint(theLine, theVertex, theArrPeriods);
120       return IntPatch_SPntSeamU;
121
122     default:
123       break;
124     }
125
126     if(anAddedPType != IntPatch_SPntNone)
127     {
128       theSingularSurfaceID = i + 1;
129       AddPointIntoLine(theLine, theArrPeriods, anApexPoint, &theVertex);
130       return anAddedPType;
131     }
132   }
133
134   return IntPatch_SPntNone;
135 }
136
137 //=======================================================================
138 //function : IntPatch_ALineToWLine
139 //purpose  : 
140 //=======================================================================
141 IntPatch_ALineToWLine::IntPatch_ALineToWLine(const Handle(Adaptor3d_HSurface)& theS1,
142                                              const Handle(Adaptor3d_HSurface)& theS2,
143                                              const Standard_Integer theNbPoints) :
144   myS1(theS1),
145   myS2(theS2),
146   myNbPointsInWline(theNbPoints),
147   myTolOpenDomain(1.e-9),
148   myTolTransition(1.e-8),
149   myTol3D(Precision::Confusion())
150
151   const GeomAbs_SurfaceType aTyps1 = theS1->GetType();
152   const GeomAbs_SurfaceType aTyps2 = theS2->GetType();
153
154   switch(aTyps1)
155   {
156   case GeomAbs_Plane:
157     myQuad1.SetValue(theS1->Plane());
158     break;
159
160   case GeomAbs_Cylinder:
161     myQuad1.SetValue(theS1->Cylinder());
162     break;
163
164   case GeomAbs_Sphere:
165     myQuad1.SetValue(theS1->Sphere());
166     break;
167
168   case GeomAbs_Cone:
169     myQuad1.SetValue(theS1->Cone());
170     break;
171
172   case GeomAbs_Torus:
173     myQuad1.SetValue(theS1->Torus());
174     break;
175
176   default:
177     break;
178   }
179
180   switch(aTyps2)
181   {
182   case GeomAbs_Plane:
183     myQuad2.SetValue(theS2->Plane());
184     break;
185   case GeomAbs_Cylinder:
186     myQuad2.SetValue(theS2->Cylinder());
187     break;
188
189   case GeomAbs_Sphere:
190     myQuad2.SetValue(theS2->Sphere());
191     break;
192
193   case GeomAbs_Cone:
194     myQuad2.SetValue(theS2->Cone());
195     break;
196
197   case GeomAbs_Torus:
198     myQuad2.SetValue(theS2->Torus());
199     break;
200
201   default:
202     break;
203   }
204 }
205
206 //=======================================================================
207 //function : SetTol3D
208 //purpose  : 
209 //=======================================================================
210 void IntPatch_ALineToWLine::SetTol3D(const Standard_Real aTol)
211 {
212   myTol3D = aTol;
213 }
214 //=======================================================================
215 //function : Tol3D
216 //purpose  : 
217 //=======================================================================
218 Standard_Real IntPatch_ALineToWLine::Tol3D()const
219 {
220   return myTol3D;
221 }
222 //=======================================================================
223 //function : SetTolTransition
224 //purpose  : 
225 //=======================================================================
226 void IntPatch_ALineToWLine::SetTolTransition(const Standard_Real aTol)
227 {
228   myTolTransition = aTol;
229 }
230 //=======================================================================
231 //function : TolTransition
232 //purpose  : 
233 //=======================================================================
234 Standard_Real IntPatch_ALineToWLine::TolTransition()const
235 {
236   return myTolTransition;
237 }
238 //=======================================================================
239 //function : SetTolOpenDomain
240 //purpose  : 
241 //=======================================================================
242 void IntPatch_ALineToWLine::SetTolOpenDomain(const Standard_Real aTol)
243 {
244   myTolOpenDomain = aTol;
245 }
246 //=======================================================================
247 //function : TolOpenDomain
248 //purpose  : 
249 //=======================================================================
250   Standard_Real IntPatch_ALineToWLine::TolOpenDomain()const
251 {
252   return myTolOpenDomain;
253 }
254
255 //=======================================================================
256 //function : GetSectionRadius
257 //purpose  : 
258 //=======================================================================
259 Standard_Real IntPatch_ALineToWLine::GetSectionRadius(const gp_Pnt& thePnt3d) const
260 {
261   Standard_Real aRetVal = RealLast();
262   for (Standard_Integer i = 0; i < 2; i++)
263   {
264     const IntSurf_Quadric& aQuad = i ? myQuad2 : myQuad1;
265     if (aQuad.TypeQuadric() == GeomAbs_Cone)
266     {
267       const gp_Cone aCone = aQuad.Cone();
268       const gp_XYZ aRVec = thePnt3d.XYZ() - aCone.Apex().XYZ();
269       const gp_XYZ &aDir = aCone.Axis().Direction().XYZ();
270
271       aRetVal = Min(aRetVal, Abs(aRVec.Dot(aDir)*Tan(aCone.SemiAngle())));
272     }
273     else if (aQuad.TypeQuadric() == GeomAbs_Sphere)
274     {
275       const gp_Sphere aSphere = aQuad.Sphere();
276       const gp_XYZ aRVec = thePnt3d.XYZ() - aSphere.Location().XYZ();
277       const gp_XYZ &aDir = aSphere.Position().Direction().XYZ();
278       const Standard_Real aR = aSphere.Radius();
279       const Standard_Real aD = aRVec.Dot(aDir);
280       const Standard_Real aDelta = aR*aR - aD*aD;
281       if (aDelta <= 0.0)
282       {
283         aRetVal = 0.0;
284         break;
285       }
286       else
287       {
288         aRetVal = Min(aRetVal, Sqrt(aDelta));
289       }
290     }
291   }
292
293   return aRetVal;
294 }
295
296 //=======================================================================
297 //function : MakeWLine
298 //purpose  : 
299 //=======================================================================
300 void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theAline,
301                                       IntPatch_SequenceOfLine& theLines) const
302
303   Standard_Boolean included;
304   Standard_Real f = theAline->FirstParameter(included); 
305   if(!included) {
306     f+=myTolOpenDomain;
307   }
308   Standard_Real l = theAline->LastParameter(included); 
309   if(!included) { 
310     l-=myTolOpenDomain;
311   }
312
313   MakeWLine(theAline, f, l, theLines);
314 }
315
316 //=======================================================================
317 //function : MakeWLine
318 //purpose  : 
319 //=======================================================================
320 void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theALine,
321                                       const Standard_Real theFPar,
322                                       const Standard_Real theLPar,
323                                       IntPatch_SequenceOfLine& theLines) const 
324 {
325   const Standard_Integer aNbVert = theALine->NbVertex();
326   if (aNbVert == 0)
327   {
328     return;
329   }
330
331 #if 0
332   //To draw ALine as a wire DRAW-object use the following code.
333   {
334     static int zzz = 0;
335     zzz++;
336
337     bool flShow = /*(zzz == 1)*/false;
338
339     if (flShow)
340     {
341       std::cout << " +++ DUMP ALine (begin) +++++" << std::endl;
342       Standard_Integer aI = 0;
343       const Standard_Real aStep = (theLPar - theFPar) / 9999.0;
344       for (Standard_Real aPrm = theFPar; aPrm < theLPar; aPrm += aStep)
345       {
346         const gp_Pnt aPP(theALine->Value(aPrm));
347         std::cout << "vertex v" << ++aI << " " << aPP.X() << " " << aPP.Y() << " " << aPP.Z() << std::endl;
348       }
349
350       gp_Pnt aPP(theALine->Value(theLPar));
351       std::cout << "vertex v" << ++aI << " " << aPP.X() << " " << aPP.Y() << " " << aPP.Z() << std::endl;
352       std::cout << " --- DUMP ALine (end) -----" << std::endl;
353     }
354   }
355
356   //Copy all output information and apply it as a TCL-code in DRAW.
357
358   //After that, use TCL-script below:
359
360   /* ********************************* Script (begin)
361   shape ww w
362   copy v1 vprev
363   for {set i 2} {$i <= 10000} {incr i} {
364     distmini dd vprev v$i;
365
366     if { [dval dd_val] > 1.0e-7} {
367       edge ee vprev v$i;
368       add ee ww;
369       copy v$i vprev;
370     }
371   }
372   ********************************** Script (end) */
373 #endif
374
375   //The same points can be marked by different vertices.
376   //The code below unifies tolerances of all vertices
377   //marking the same points.
378   for (Standard_Integer i = 1; i < aNbVert; i++)
379   {
380     IntPatch_Point &aCurVert = theALine->ChangeVertex(i);
381     const IntSurf_PntOn2S &aCurrPt = aCurVert.PntOn2S();
382     const Standard_Real aCurToler = aCurVert.Tolerance();
383     for (Standard_Integer j = i + 1; j <= aNbVert; j++)
384     {
385       IntPatch_Point &aVert = theALine->ChangeVertex(j);
386       const IntSurf_PntOn2S &aNewPt = aVert.PntOn2S();
387       const Standard_Real aToler = aVert.Tolerance();
388
389       const Standard_Real aSumTol = aCurToler + aToler;
390       if (aCurrPt.IsSame(aNewPt, aSumTol))
391       {
392         aCurVert.SetTolerance(aSumTol);
393         aVert.SetTolerance(aSumTol);
394       }
395     }
396   }
397
398   const Standard_Real aTol = 2.0*myTol3D+Precision::Confusion();
399   const Standard_Real aPrmTol = Max(1.0e-4*(theLPar - theFPar), Precision::PConfusion());
400
401   IntPatch_SpecPntType aPrePointExist = IntPatch_SPntNone;
402   
403   NCollection_Array1<Standard_Real> aVertexParams(1, aNbVert);
404   NCollection_Array1<IntPatch_Point> aSeqVertex(1, aNbVert);
405
406   //It is possible to have several vertices with equal parameters.
407   NCollection_Array1<Standard_Boolean> hasVertexBeenChecked(1, aNbVert);
408
409   Handle(IntSurf_LineOn2S) aLinOn2S;
410   Standard_Real aParameter = theFPar;
411   
412   for(Standard_Integer i = aVertexParams.Lower(); i <= aVertexParams.Upper(); i++)
413   {
414     const IntPatch_Point& aVert = theALine->Vertex(i);
415     const Standard_Real aPar = aVert.ParameterOnLine();
416     aVertexParams(i) = aPar;
417     hasVertexBeenChecked(i) = Standard_False;
418   }
419
420   Standard_Integer aSingularSurfaceID = 0;
421   Standard_Real anArrPeriods[] = { 0.0,  //U1
422                                    0.0,  //V1
423                                    0.0,  //U2
424                                    0.0}; //V2
425
426   IntSurf::SetPeriod(myS1, myS2, anArrPeriods);
427
428   IntSurf_PntOn2S aPrevLPoint;
429
430   while(aParameter < theLPar)
431   {
432     Standard_Real aStep = (theLPar - aParameter) / (Standard_Real)(myNbPointsInWline - 1);
433     if(aStep < Epsilon(theLPar))
434       break;
435
436     Standard_Integer aNewVertID = 0;
437     aLinOn2S = new IntSurf_LineOn2S;
438     
439     const Standard_Real aStepMin = 0.1*aStep, aStepMax = 10.0*aStep;
440
441     Standard_Boolean isLast = Standard_False;
442     Standard_Real aPrevParam = aParameter;
443     for(; !isLast; aParameter += aStep)
444     {
445       IntSurf_PntOn2S aPOn2S;
446
447       if(theLPar <= aParameter)
448       {
449         isLast = Standard_True;
450         if(aPrePointExist != IntPatch_SPntNone)
451         {
452           break;
453         }
454         else
455         {
456           aParameter = theLPar;
457         }
458       }
459
460       Standard_Boolean isPointValid = Standard_False;
461       Standard_Real aTgMagn = 0.0;
462       {
463         gp_Pnt aPnt3d;
464         gp_Vec aTg;
465         theALine->D1(aParameter, aPnt3d, aTg);
466         if (GetSectionRadius(aPnt3d) < 5.0e-6)
467         {
468           // We cannot compute 2D-parameters of
469           // aPOn2S correctly.
470
471           isPointValid = Standard_False;
472         }
473         else
474         {
475           isPointValid = Standard_True;
476         }
477         
478         aTgMagn = aTg.Magnitude();
479         Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;
480         myQuad1.Parameters(aPnt3d, u1, v1);
481         myQuad2.Parameters(aPnt3d, u2, v2);
482         aPOn2S.SetValue(aPnt3d, u1, v1, u2, v2);
483       }
484
485       if(aPrePointExist != IntPatch_SPntNone)
486       {
487         const Standard_Real aURes = Max(myS1->UResolution(myTol3D),
488                                                 myS2->UResolution(myTol3D)),
489                             aVRes = Max(myS1->VResolution(myTol3D),
490                                                 myS2->VResolution(myTol3D));
491
492         const Standard_Real aTol2d = (aPrePointExist == IntPatch_SPntPole) ? -1.0 : 
493                     (aPrePointExist == IntPatch_SPntSeamV)? aVRes :
494                     (aPrePointExist == IntPatch_SPntSeamUV)? Max(aURes, aVRes) : aURes;
495
496         IntSurf_PntOn2S aRPT = aPOn2S;
497
498         if (aPrePointExist == IntPatch_SPntPole)
499         {
500           Standard_Real aPrt = 0.5*(aPrevParam + theLPar);
501           for (Standard_Integer i = aVertexParams.Lower();
502                i <= aVertexParams.Upper(); i++)
503           {
504             const Standard_Real aParam = aVertexParams(i);
505
506             if (aParam <= aPrevParam)
507               continue;
508
509             if ((aParam - aPrevParam) < aPrmTol)
510             {
511               const gp_Pnt aPnt3d(theALine->Value(aParam));
512               if (aPOn2S.Value().SquareDistance(aPnt3d) < Precision::SquareConfusion())
513               {
514                 // i-th vertex is the same as a Pole/Apex.
515                 // So, it should be ignored.
516                 continue;
517               }
518             }
519
520             aPrt = 0.5*(aParam + aPrevParam);
521             break;
522           }
523
524           const gp_Pnt aPnt3d(theALine->Value(aPrt));
525           Standard_Real u1, v1, u2, v2;
526           myQuad1.Parameters(aPnt3d, u1, v1);
527           myQuad2.Parameters(aPnt3d, u2, v2);
528           aRPT.SetValue(aPnt3d, u1, v1, u2, v2);
529
530           if (aPOn2S.IsSame(aPrevLPoint, Max(Precision::Approximation(), aTol)))
531           {
532             //Set V-parameter as precise value found on the previous step.
533             if (aSingularSurfaceID == 1)
534             {
535               aPOn2S.ParametersOnS1(u2, v2);
536               aPOn2S.SetValue(Standard_True, u1, v2);
537             }
538             else //if (aSingularSurfaceID == 2)
539             {
540               aPOn2S.ParametersOnS2(u1, v1);
541               aPOn2S.SetValue(Standard_False, u2, v1);
542             }
543           }
544         }
545
546         if(IntPatch_SpecialPoints::
547                       ContinueAfterSpecialPoint(myS1, myS2, aRPT,
548                                                 aPrePointExist, aTol2d,
549                                                 aPrevLPoint, Standard_False))
550         {
551           AddPointIntoLine(aLinOn2S, anArrPeriods, aPrevLPoint);
552         }
553         else if(aParameter == theLPar)
554         {// Strictly equal!!!
555           break;
556         }
557       }
558
559       aPrePointExist = IntPatch_SPntNone;
560
561       Standard_Integer aVertexNumber = -1;
562       for(Standard_Integer i = aVertexParams.Lower(); i <= aVertexParams.Upper(); i++)
563       {
564         if(hasVertexBeenChecked(i))
565           continue;
566
567         const IntPatch_Point &aVP = theALine->Vertex(i);
568         const Standard_Real aParam = aVertexParams(i);
569         if( ((aPrevParam < aParam) && (aParam <= aParameter)) ||
570             ((aPrevParam == aParameter) && (aParam == aParameter))||
571             (aPOn2S.IsSame(aVP.PntOn2S(), aVP.Tolerance()) && 
572                     (Abs(aVP.ParameterOnLine() - aParameter) < aPrmTol)))
573         {
574           //We have either jumped over the vertex or "fell" on the vertex.
575           //However, ALine can be self-interfered. Therefore, we need to check
576           //vertex parameter and 3D-distance together.
577
578           aVertexNumber = i;
579           break;
580         }
581       }
582
583       aPrevParam = aParameter;
584       
585       if(aVertexNumber < 0)
586       {
587         if (isPointValid)
588         {
589           StepComputing(theALine, aPOn2S, theLPar, aParameter, aTgMagn,
590                         aStepMin, aStepMax, myTol3D, aStep);
591           AddPointIntoLine(aLinOn2S, anArrPeriods, aPOn2S);
592           aPrevLPoint = aPOn2S;
593         }
594
595         continue;
596       }
597
598       IntPatch_Point aVtx = theALine->Vertex(aVertexNumber);
599       const Standard_Real aNewVertexParam = aLinOn2S->NbPoints() + 1;
600
601       //ATTENTION!!!
602       // IsPoleOrSeam inserts new point in aLinOn2S if aVtx respects
603       //to some special point. Otherwise, aLinOn2S is not changed.
604
605       // Find a point for reference parameter. It will be used
606       // if real parameter value cannot be precise (see comment to 
607       // IsPoleOrSeam(...) function). 
608       IntSurf_PntOn2S aPrefIso = aVtx.PntOn2S();
609       if (aLinOn2S->NbPoints() < 1)
610       {
611         for (Standard_Integer i = aVertexNumber + 1; i <= aVertexParams.Upper(); i++)
612         {
613           const Standard_Real aParam = aVertexParams(i);
614           if ((aParam - aVertexParams(aVertexNumber)) > Precision::PConfusion())
615           {
616             const Standard_Real aPrm = 0.5*(aParam + aVertexParams(aVertexNumber));
617             const gp_Pnt aPnt3d(theALine->Value(aPrm));
618             Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;
619             myQuad1.Parameters(aPnt3d, u1, v1);
620             myQuad2.Parameters(aPnt3d, u2, v2);
621             aPrefIso.SetValue(aPnt3d, u1, v1, u2, v2);
622             break;
623           }
624         }
625       }
626       else
627       {
628         aPrefIso = aLinOn2S->Value(aLinOn2S->NbPoints());
629       }
630
631       aPrePointExist = IsPoleOrSeam(myS1, myS2, aPrefIso, aLinOn2S, aVtx,
632                                 anArrPeriods, aTol, aSingularSurfaceID);
633
634       const Standard_Real aCurVertParam = aVtx.ParameterOnLine();
635       if(aPrePointExist != IntPatch_SPntNone)
636       {
637         aPrevParam = aParameter = aCurVertParam;
638       }
639       else
640       {
641         if (!isPointValid)
642         {
643           //Take a farther point of ALine (with greater parameter)
644           continue;
645         }
646
647         if(aVtx.Tolerance() > aTol)
648         {
649           aVtx.SetValue(aPOn2S);
650           AddPointIntoLine(aLinOn2S, anArrPeriods, aPOn2S);
651         }
652         else
653         {
654           AddVertexPoint(aLinOn2S, aVtx, anArrPeriods);
655         }
656       }
657
658       aPrevLPoint = aPOn2S = aLinOn2S->Value(aLinOn2S->NbPoints());
659
660       {
661         Standard_Boolean isFound = Standard_False;
662         const Standard_Real aSqTol = aTol*aTol;
663         const gp_Pnt aP1(theALine->Value(aCurVertParam));
664         const IntSurf_PntOn2S& aVertP2S = aVtx.PntOn2S();
665         const Standard_Real aVertToler  = aVtx.Tolerance();
666
667         for(Standard_Integer i = aVertexParams.Lower(); i <= aVertexParams.Upper(); i++)
668         {
669           if(hasVertexBeenChecked(i))
670             continue;
671
672           const gp_Pnt aP2(theALine->Value(aVertexParams(i)));
673           
674           if(aP1.SquareDistance(aP2) < aSqTol)
675           {
676             IntPatch_Point aLVtx = theALine->Vertex(i);
677             aLVtx.SetValue(aVertP2S);
678             aLVtx.SetTolerance(aVertToler);
679             aLVtx.SetParameter(aNewVertexParam);
680             aSeqVertex(++aNewVertID) = aLVtx;
681             hasVertexBeenChecked(i) = Standard_True;
682             isFound = Standard_True;
683           }
684           else if(isFound)
685           {
686             break;
687           }
688         }
689       }
690
691       if ((aPrePointExist != IntPatch_SPntNone) && (aLinOn2S->NbPoints() > 1))
692         break;
693     }//for(; !isLast; aParameter += aStep)
694
695     if(aLinOn2S->NbPoints() < 2)
696     {
697       aParameter += aStep;
698       continue;
699     }
700
701     //-----------------------------------------------------------------
702     //--  Computation of transitions of the line on two surfaces    ---
703     //-----------------------------------------------------------------
704     IntSurf_TypeTrans trans1,trans2;
705     {
706       Standard_Integer indice1;
707       Standard_Real dotcross;
708       gp_Pnt aPP0, aPP1;
709       //
710       trans1=IntSurf_Undecided;
711       trans2=IntSurf_Undecided;
712       //
713       indice1 = aLinOn2S->NbPoints()/3;
714       if(indice1<=2) {
715         indice1 = 2;
716       }
717       //
718       aPP1=aLinOn2S->Value(indice1).Value();
719       aPP0=aLinOn2S->Value(indice1-1).Value();
720       //
721       gp_Vec tgvalid(aPP0, aPP1);
722       gp_Vec aNQ1 = myQuad1.Normale(aPP0);
723       gp_Vec aNQ2 = myQuad2.Normale(aPP0);
724       //
725       dotcross = tgvalid.DotCross(aNQ2, aNQ1);
726       if (dotcross > myTolTransition) {
727         trans1 = IntSurf_Out;
728         trans2 = IntSurf_In;
729       }
730       else if(dotcross < -myTolTransition) { 
731         trans1 = IntSurf_In;
732         trans2 = IntSurf_Out;
733       } 
734     }
735
736     //-----------------------------------------------------------------
737     //--              W  L  i  n  e     c  r  e  a  t  i  o  n      ---
738     //-----------------------------------------------------------------
739     Handle(IntPatch_WLine) aWLine;
740     //
741     if(theALine->TransitionOnS1() == IntSurf_Touch)  { 
742                                   aWLine = new IntPatch_WLine(aLinOn2S,
743                                   theALine->IsTangent(),
744                                   theALine->SituationS1(),
745                                   theALine->SituationS2());
746     }
747     else if(theALine->TransitionOnS1() == IntSurf_Undecided)  { 
748       aWLine = new IntPatch_WLine(aLinOn2S, theALine->IsTangent());
749     }
750     else { 
751       aWLine = new IntPatch_WLine(aLinOn2S, theALine->IsTangent(),
752                                   trans1, // aline->TransitionOnS1(),
753                                   trans2);  //aline->TransitionOnS2());
754     }
755
756     for(Standard_Integer i = aSeqVertex.Lower(); i <= aNewVertID; i++)
757     {
758       const IntPatch_Point& aVtx = aSeqVertex(i);
759       aWLine->AddVertex(aVtx);
760     }
761
762     aWLine->SetPeriod(anArrPeriods[0],anArrPeriods[1],anArrPeriods[2],anArrPeriods[3]);
763
764     //the method ComputeVertexParameters can reduce the number of points in <aWLine>
765     aWLine->ComputeVertexParameters(myTol3D);
766
767     if (aWLine->NbPnts() > 1)
768     {
769       aWLine->EnablePurging(Standard_False);
770 #ifdef INTPATCH_ALINETOWLINE_DEBUG
771       aWLine->Dump(0);
772 #endif
773       theLines.Append(aWLine);
774     }
775   }//while(aParameter < theLPar)
776 }
777
778 //=======================================================================
779 //function : CheckDeflection
780 //purpose  : Returns:
781 //            -1 - step is too small
782 //            0  - step is normal
783 //            +1 - step is too big
784 //=======================================================================
785 Standard_Integer IntPatch_ALineToWLine::CheckDeflection(const gp_XYZ& theMidPt,
786                                                         const Standard_Real theMaxDeflection) const
787 {
788   Standard_Real aDist = Abs(myQuad1.Distance(theMidPt));
789   if(aDist > theMaxDeflection)
790     return 1;
791
792   aDist = Max(Abs(myQuad2.Distance(theMidPt)), aDist);
793   
794   if(aDist > theMaxDeflection)
795     return 1;
796
797   if((aDist + aDist) < theMaxDeflection)
798     return -1;
799
800   return 0;
801 }
802
803 //=======================================================================
804 //function : StepComputing
805 //purpose  : 
806 //=======================================================================
807 Standard_Boolean IntPatch_ALineToWLine::
808                       StepComputing(const Handle(IntPatch_ALine)& theALine,
809                                     const IntSurf_PntOn2S& thePOn2S,
810                                     const Standard_Real theLastParOfAline,
811                                     const Standard_Real theCurParam,
812                                     const Standard_Real theTgMagnitude,
813                                     const Standard_Real theStepMin,
814                                     const Standard_Real theStepMax,
815                                     const Standard_Real theMaxDeflection,
816                                     Standard_Real& theStep) const
817 {
818   if(theTgMagnitude < Precision::Confusion())
819     return Standard_False;
820   
821   const Standard_Real anEps = myTol3D;
822
823   //Indeed, 1.0e+15 < 2^50 < 1.0e+16. Therefore,
824   //if we apply bisection method to the range with length
825   //1.0e+6 then we will be able to find solution with max error ~1.0e-9.
826   const Standard_Integer aNbIterMax = 50;
827
828   const Standard_Real aNotFilledRange = theLastParOfAline - theCurParam;
829   Standard_Real aMinStep = theStepMin, aMaxStep = Min(theStepMax, aNotFilledRange);
830
831   if(aMinStep > aMaxStep)
832   {
833     theStep = aMaxStep;
834     return Standard_True;
835   }
836
837   const Standard_Real aR = IntPatch_PointLine::
838                             CurvatureRadiusOfIntersLine(myS1, myS2, thePOn2S);
839
840 #if 0
841   {
842     static int zzz = 0;
843     zzz++;
844     std::cout << "*** R" << zzz << " (begin)" << std::endl;
845     Standard_Real aU1, aV1, aU2, aV2;
846     thePOn2S.Parameters(aU1, aV1, aU2, aV2);
847     std::cout << "Prms: " << aU1 << ", " << aV1 << ", " << aU2 << ", " << aV2 << std::endl;
848     std::cout << "Radius = " << aR << std::endl;
849     std::cout << "*** R" << zzz << " (end)" << std::endl;
850   }
851 #endif
852
853   if(aR < 0.0)
854   {
855     return Standard_False;
856   }
857   else
858   {
859     //The 3D-step is defined as length of the tangent to the osculating circle
860     //by the condition that the distance from end point of the tangent to the
861     //circle is no greater than anEps. theStep is the step in
862     //parameter space of intersection curve (must be converted from 3D-step).
863
864     theStep = Min(sqrt(anEps*(2.0*aR + anEps))/theTgMagnitude, aMaxStep);
865     theStep = Max(theStep, aMinStep);
866   }
867
868   //The step value has been computed for osculating circle.
869   //Now it should be checked for real intersection curve
870   //and is made more precise in case of necessity.
871
872   Standard_Integer aNbIter = 0;
873   do
874   {
875     aNbIter++;
876
877     const gp_XYZ& aP1 = thePOn2S.Value().XYZ();
878     const gp_XYZ aP2(theALine->Value(theCurParam + theStep).XYZ());
879     const Standard_Integer aStatus = CheckDeflection(0.5*(aP1 + aP2), theMaxDeflection);
880
881     if(aStatus == 0)
882       break;
883
884     if(aStatus < 0)
885     {
886       aMinStep = theStep;
887     }
888     else //if(aStatus > 0)
889     {
890       aMaxStep = theStep;
891     }
892
893     theStep = 0.5*(aMinStep + aMaxStep);
894   }
895   while(((aMaxStep - aMinStep) > Precision::PConfusion()) && (aNbIter <= aNbIterMax));
896
897   if(aNbIter > aNbIterMax)
898     return Standard_False;
899
900   return Standard_True;
901 }