0031007: Coding - eliminate warnings issued while compiling with -pedantic flag
[occt.git] / src / Bisector / Bisector_BisecCC.cxx
1 // Created on: 1994-03-10
2 // Created by: Yves FRICAUD
3 // Copyright (c) 1994-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
18 #include <Bisector.hxx>
19 #include <Bisector_BisecCC.hxx>
20 #include <Bisector_BisecPC.hxx>
21 #include <Bisector_Curve.hxx>
22 #include <Bisector_FunctionH.hxx>
23 #include <Bisector_PointOnBis.hxx>
24 #include <Bisector_PolyBis.hxx>
25 #include <Geom2d_Circle.hxx>
26 #include <Geom2d_Curve.hxx>
27 #include <Geom2d_Geometry.hxx>
28 #include <Geom2d_Line.hxx>
29 #include <Geom2d_TrimmedCurve.hxx>
30 #include <Geom2dAdaptor_Curve.hxx>
31 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
32 #include <Geom2dGcc.hxx>
33 #include <Geom2dGcc_Circ2d2TanRad.hxx>
34 #include <Geom2dGcc_QualifiedCurve.hxx>
35 #include <Geom2dInt_GInter.hxx>
36 #include <Geom2dLProp_CLProps2d.hxx>
37 #include <gp.hxx>
38 #include <gp_Pnt2d.hxx>
39 #include <gp_Trsf2d.hxx>
40 #include <gp_Vec2d.hxx>
41 #include <IntRes2d_IntersectionPoint.hxx>
42 #include <math_BissecNewton.hxx>
43 #include <math_FunctionRoot.hxx>
44 #include <math_FunctionRoots.hxx>
45 #include <Precision.hxx>
46 #include <Standard_DivideByZero.hxx>
47 #include <Standard_DomainError.hxx>
48 #include <Standard_NotImplemented.hxx>
49 #include <Standard_OutOfRange.hxx>
50 #include <Standard_RangeError.hxx>
51 #include <Standard_Type.hxx>
52
53 IMPLEMENT_STANDARD_RTTIEXT(Bisector_BisecCC,Bisector_Curve)
54
55 static Standard_Boolean ProjOnCurve (const gp_Pnt2d& P,
56   const Handle(Geom2d_Curve)& C,
57   Standard_Real& theParam);
58
59 static Standard_Real    Curvature (const Handle(Geom2d_Curve)& C,
60   Standard_Real         U,
61   Standard_Real         Tol) ;
62
63 static Standard_Boolean TestExtension (const Handle(Geom2d_Curve)& C1,
64   const Handle(Geom2d_Curve)& C2,
65   const Standard_Integer      Start_End);
66
67 static Standard_Boolean DiscretPar(const Standard_Real     DU, 
68   const Standard_Real     EpsMin, 
69   const Standard_Real     EpsMax,
70   const Standard_Integer  NbMin,
71   const Standard_Integer  NbMax,
72   Standard_Real&    Eps,
73   Standard_Integer& Nb);
74
75 //=============================================================================
76 //function :
77 //purpose  :
78 //=============================================================================
79 Bisector_BisecCC::Bisector_BisecCC()
80 {
81   shiftParameter = 0;
82   isEmpty        = Standard_False;
83 }
84
85 //=============================================================================
86 //function :
87 //purpose  :
88 //=============================================================================
89 Bisector_BisecCC::Bisector_BisecCC(const Handle(Geom2d_Curve)& Cu1,
90   const Handle(Geom2d_Curve)& Cu2,
91   const Standard_Real         Side1,
92   const Standard_Real         Side2,
93   const gp_Pnt2d&             Origin,
94   const Standard_Real         DistMax)
95 {
96   Perform (Cu1,Cu2,Side1,Side2,Origin,DistMax);
97 }
98
99 //=============================================================================
100 //function : Perform
101 //purpose  :
102 //=============================================================================
103 void  Bisector_BisecCC::Perform(const Handle(Geom2d_Curve)& Cu1,
104   const Handle(Geom2d_Curve)& Cu2,
105   const Standard_Real         Side1,
106   const Standard_Real         Side2,
107   const gp_Pnt2d&             Origin,                           
108   const Standard_Real         DistMax)
109 {
110   isEmpty = Standard_False;
111   distMax = DistMax;
112
113   curve1  = Handle (Geom2d_Curve)::DownCast(Cu1->Copy());
114   curve2  = Handle (Geom2d_Curve)::DownCast(Cu2->Copy());
115
116   sign1 = Side1;
117   sign2 = Side2;
118
119   isConvex1 = Bisector::IsConvex(curve1,sign1);
120   isConvex2 = Bisector::IsConvex(curve2,sign2);
121
122   Standard_Real    U,UC1,UC2,Dist,dU,USol;  
123   gp_Pnt2d         P;
124   Standard_Integer NbPnts    = 21;
125   Standard_Real    EpsMin    = 10*Precision::Confusion();
126   Standard_Boolean YaPoly    = Standard_True;
127   Standard_Boolean OriInPoly = Standard_False;
128   //---------------------------------------------
129   // Calculate first point of the polygon.
130   //---------------------------------------------
131   Standard_Boolean isProjDone = ProjOnCurve (Origin,curve1, U);
132
133   if(!isProjDone)
134   {
135     isEmpty = Standard_True;
136     return;
137   }
138
139   P = ValueByInt  (U,UC1,UC2,Dist);
140   if(Dist < Precision::Confusion())
141   {
142     gp_Pnt2d aP1 = curve1->Value(UC1);
143     gp_Pnt2d aP2 = curve2->Value(UC2);
144     Standard_Real dp = (aP1.Distance(P)+aP2.Distance(P));
145     Standard_Real dorig = (aP1.Distance(Origin)+aP2.Distance(Origin));
146     if(dp < dorig)
147     {
148       isEmpty = Standard_True;
149       return;
150     }
151   }
152
153   if (Dist < Precision::Infinite()) {
154     //----------------------------------------------------
155     // the parameter of the origin point gives a point 
156     // on the polygon.
157     //----------------------------------------------------
158     myPolygon.Append(Bisector_PointOnBis(UC1,UC2,U,Dist,P));     
159     startIntervals.Append(U);
160     if (P.IsEqual(Origin,Precision::Confusion())) {    
161       //----------------------------------------
162       // test if the first point is the origin.
163       //----------------------------------------
164       OriInPoly = Standard_True;
165     }
166   }
167   else {
168     //-------------------------------------------------------
169     // The origin point is on the extension.
170     // Find the first point of the polygon by dichotomy.
171     //-------------------------------------------------------
172     dU     = (curve1->LastParameter() - U)/(NbPnts - 1); 
173     U += dU;
174     for (Standard_Integer i = 1; i <= NbPnts - 1; i++) {
175       P    = ValueByInt(U,UC1,UC2,Dist);
176       if (Dist < Precision::Infinite()) {
177         USol = SearchBound(U - dU,U);    
178         P    = ValueByInt(USol,UC1,UC2,Dist);   
179         startIntervals.Append(USol);
180         myPolygon.Append(Bisector_PointOnBis(UC1,UC2,USol,Dist,P));
181         break;
182       }
183       U += dU;
184     }
185   }
186
187   if ( myPolygon.Length() != 0 ) {
188     SupLastParameter();
189     //----------------------------------------------
190     // Construction of the polygon of the bissectrice.
191     //---------------------------------------------
192     U                = FirstParameter();
193     Standard_Real DU = LastParameter() - U;
194
195     if (DU < EpsMin) {NbPnts = 3;} 
196     dU = DU/(NbPnts - 1);
197
198     U += dU;
199     //  modified by NIZHNY-EAP Fri Jan 21 09:33:20 2000 ___BEGIN___
200     //  prevent addition of the same point
201     gp_Pnt2d prevPnt = P;
202     for (Standard_Integer i = 1; i <= NbPnts - 1; i++) {
203       P    = ValueByInt(U,UC1,UC2,Dist);
204       if (Dist < Precision::Infinite()) {
205         if (P.Distance (prevPnt) > Precision::Confusion())
206           myPolygon.Append(Bisector_PointOnBis(UC1,UC2,U,Dist,P));
207       }
208       else {
209         USol = SearchBound(U - dU,U);
210         P    = ValueByInt(USol,UC1,UC2,Dist);      
211         endIntervals.SetValue(1,USol);
212         if (P.Distance (prevPnt) > Precision::Confusion())
213           myPolygon.Append(Bisector_PointOnBis(UC1,UC2,USol,Dist,P));
214         break;
215       }
216       U += dU;
217       prevPnt=P;
218       //  modified by NIZHNY-EAP Fri Jan 21 09:33:24 2000 ___END___
219     }  
220   }
221   else {
222     //----------------
223     // Empty Polygon.    
224     //----------------
225     YaPoly = Standard_False;
226   }
227
228   extensionStart = Standard_False;
229   extensionEnd   = Standard_False;
230   pointStart     = Origin;
231
232   if (isConvex1 && isConvex2) {
233     if (YaPoly) pointEnd = myPolygon.Last().Point();
234   }
235   else {
236     //-----------------------------------------------------------------------------
237     // Extension : The curve is extended at the beginning and/or the end if
238     //                - one of two curves is concave.
239     //                - the curves have a common point at the beginning and/or the end
240     //                - the angle of opening at the common point between two curves
241     //                  values M_PI.
242     // the extension at the beginning  is taken into account if the origin is found above.
243     // ie : the origin is not the in the polygon.
244     //-----------------------------------------------------------------------------
245
246     //---------------------------------
247     // Do the extensions exist ?
248     //---------------------------------
249     if (OriInPoly) {
250       extensionStart = Standard_False;
251     }
252     else {
253       extensionStart = TestExtension(curve1,curve2,1);
254     }
255     extensionEnd     = TestExtension(curve1,curve2,2);
256
257     //-----------------
258     // Calculate pointEnd.
259     //-----------------
260     if (extensionEnd) {
261       pointEnd = curve1->Value(curve1->LastParameter());
262     }
263     else if (YaPoly) {
264       pointEnd = myPolygon.Last().Point();
265     }
266     else {
267       ComputePointEnd();
268     }
269     //------------------------------------------------------
270     // Update the Limits of intervals of  definition.
271     //------------------------------------------------------
272     if (YaPoly) {
273       if (extensionStart) {
274         gp_Pnt2d       P1     = myPolygon.First().Point();
275         Standard_Real  UFirst = startIntervals.First() - pointStart.Distance(P1);
276         startIntervals.InsertBefore(1,UFirst);
277         endIntervals  .InsertBefore(1,startIntervals.Value(2));
278       }
279       if (extensionEnd) {
280         gp_Pnt2d       P1;
281         Standard_Real  UFirst,ULast;
282         P1     = myPolygon.Last().Point();
283         UFirst = endIntervals.Last();      
284         ULast  = UFirst + pointEnd.Distance(P1);
285         startIntervals.Append(UFirst);
286         endIntervals  .Append(ULast );
287       }
288     }
289     else {
290       //--------------------------------------------------
291       // No polygon => the bissectrice is a segment.
292       //--------------------------------------------------
293       startIntervals.Append(0.);
294       endIntervals  .Append(pointEnd.Distance(pointStart));
295     }
296   }
297   if (!YaPoly && !extensionStart && !extensionEnd) 
298     isEmpty = Standard_True;
299   //  modified by NIZHNY-EAP Mon Jan 17 17:32:40 2000 ___BEGIN___
300   if (myPolygon.Length() <= 2)
301     isEmpty = Standard_True;
302   //  modified by NIZHNY-EAP Mon Jan 17 17:32:42 2000 ___END___
303 }
304
305 //=============================================================================
306 //function : IsExtendAtStart
307 //purpose  :
308 //=============================================================================
309 Standard_Boolean Bisector_BisecCC::IsExtendAtStart() const
310 {
311   return extensionStart;
312 }
313
314 //=============================================================================
315 //function : IsExtendAtEnd
316 //purpose  :
317 //=============================================================================
318 Standard_Boolean Bisector_BisecCC::IsExtendAtEnd() const
319 {
320   return extensionEnd;
321 }
322
323 //=============================================================================
324 //function : IsEmpty
325 //purpose :
326 //=============================================================================
327 Standard_Boolean Bisector_BisecCC::IsEmpty() const
328 {
329   return isEmpty;
330 }
331
332 //=============================================================================
333 //function : Reverse
334 //purpose  :
335 //=============================================================================
336 void Bisector_BisecCC::Reverse()
337 {
338   throw Standard_NotImplemented();
339 }
340
341 //=============================================================================
342 //function : ReversedParameter
343 // purpose :
344 //=============================================================================
345 Standard_Real  Bisector_BisecCC::ReversedParameter(const Standard_Real U) const
346
347   return LastParameter() + FirstParameter() - U;
348 }
349
350 //=============================================================================
351 //function : Copy
352 //purpose  :
353 //=============================================================================
354 Handle(Geom2d_Geometry) Bisector_BisecCC::Copy() const
355 {
356   Handle(Geom2d_Curve) CopyCurve1
357     = Handle(Geom2d_Curve)::DownCast(curve1->Copy());
358   Handle(Geom2d_Curve) CopyCurve2
359     = Handle(Geom2d_Curve)::DownCast(curve2->Copy());
360
361   Handle(Bisector_BisecCC) C = new Bisector_BisecCC();
362
363   C -> Curve    (1, CopyCurve1)   ; C -> Curve    (2, CopyCurve2);
364   C -> Sign     (1, sign1 )       ; C -> Sign     (2, sign2 );
365   C -> IsConvex (1, isConvex1)    ; C -> IsConvex (2, isConvex2);
366   C -> Polygon  (myPolygon);
367   C -> IsEmpty  (isEmpty)  ;
368   C -> DistMax  (distMax)  ;
369   C -> StartIntervals (startIntervals); C -> EndIntervals (endIntervals);
370   C -> ExtensionStart (extensionStart); C -> ExtensionEnd (extensionEnd);
371   C -> PointStart     (pointStart)    ; C -> PointEnd     (pointEnd)    ;
372
373   return C;  
374 }
375
376 //=============================================================================
377 //function : ChangeGuide
378 //purpose  : Changet of the guideline for the parameters of the bissectrice
379 //           ATTENTION : - This can invert the direction  of parameterization.
380 //                       - This concerns only the part of the curve
381 //                         corresponding to the polygon.
382 //=============================================================================
383 Handle(Bisector_BisecCC) Bisector_BisecCC::ChangeGuide() const
384 {
385   Handle(Bisector_BisecCC) C = new Bisector_BisecCC();
386
387   C -> Curve    (1, curve2)   ; C -> Curve    (2, curve1);
388   C -> Sign     (1, sign2 )   ; C -> Sign     (2, sign1 );
389   C -> IsConvex (1, isConvex2); C -> IsConvex (2, isConvex1); 
390
391   //-------------------------------------------------------------------------
392   // Construction of the new polygon from the initial one.
393   // inversion of PointOnBis and Calculation of new parameters on the bissectrice.
394   //-------------------------------------------------------------------------
395   Bisector_PolyBis Poly;
396   if (sign1 == sign2 ) {
397     //---------------------------------------------------------------
398     // elements of the new polygon are ranked in the other direction.
399     //---------------------------------------------------------------
400     for (Standard_Integer i = myPolygon.Length(); i >=1; i--) {
401       Bisector_PointOnBis P = myPolygon.Value(i);
402       Bisector_PointOnBis NewP (P.ParamOnC2(), P.ParamOnC1(), 
403         P.ParamOnC2(), P.Distance (),
404         P.Point());
405       Poly.Append (NewP);
406     }
407   }
408   else {
409     for (Standard_Integer i = 1; i <= myPolygon.Length(); i ++) {      
410       Bisector_PointOnBis P = myPolygon.Value(i);
411       Bisector_PointOnBis NewP (P.ParamOnC2(), P.ParamOnC1(), 
412         P.ParamOnC2(), P.Distance (),
413         P.Point());
414       Poly.Append (NewP);
415     }
416   }
417   C -> Polygon        (Poly);
418   C -> FirstParameter (Poly.First().ParamOnBis());
419   C -> LastParameter  (Poly.Last() .ParamOnBis());
420
421   return C;  
422 }
423
424 //=============================================================================
425 //function : Transform
426 //purpose  :
427 //=============================================================================
428 void Bisector_BisecCC::Transform (const gp_Trsf2d& T)
429 {
430   curve1    ->Transform(T);
431   curve2    ->Transform(T);
432   myPolygon . Transform(T);
433   pointStart. Transform(T);
434   pointEnd  . Transform(T);
435 }
436
437 //=============================================================================
438 //function : IsCN
439 //purpose  :
440 //=============================================================================
441 Standard_Boolean Bisector_BisecCC::IsCN (const Standard_Integer N) const 
442 {
443   return (curve1->IsCN(N+1) && curve2->IsCN(N+1));
444 }
445
446 //=============================================================================
447 //function : FirstParameter
448 //purpose :
449 //=============================================================================
450 Standard_Real Bisector_BisecCC::FirstParameter() const
451 {
452   return startIntervals.First();
453 }
454
455 //=============================================================================
456 //function : LastParameter
457 //purpose  :
458 //=============================================================================
459 Standard_Real Bisector_BisecCC::LastParameter() const
460 {
461   return endIntervals.Last();
462 }
463
464 //=============================================================================
465 //function : Continuity
466 //purpose  :
467 //=============================================================================
468 GeomAbs_Shape Bisector_BisecCC::Continuity() const 
469 {
470   GeomAbs_Shape Cont = curve1->Continuity();
471   switch (Cont) {
472   case GeomAbs_C1 : return GeomAbs_C0; 
473   case GeomAbs_C2 : return GeomAbs_C1;
474   case GeomAbs_C3 : return GeomAbs_C2;
475   case GeomAbs_CN : return GeomAbs_CN;  
476   default: break;
477   }
478   return GeomAbs_C0;
479 }
480
481 //=============================================================================
482 //function : NbIntervals
483 //purpose  :
484 //=============================================================================
485 Standard_Integer Bisector_BisecCC::NbIntervals() const
486 {
487   return startIntervals.Length();
488 }
489
490 //=============================================================================
491 //function : IntervalFirst
492 //purpose  :
493 //=============================================================================
494 Standard_Real Bisector_BisecCC::IntervalFirst(const Standard_Integer Index) const
495 {
496   return startIntervals.Value(Index);
497 }
498
499 //=============================================================================
500 //function : IntervalLast
501 //purpose  :
502 //=============================================================================
503 Standard_Real Bisector_BisecCC::IntervalLast(const Standard_Integer Index) const
504 {
505   return endIntervals.Value(Index);
506 }
507
508 //=============================================================================
509 //function : IntervalContinuity
510 //purpose  :
511 //=============================================================================
512 GeomAbs_Shape Bisector_BisecCC::IntervalContinuity() const
513 {
514   GeomAbs_Shape Cont = curve1->Continuity();
515   switch (Cont) {
516   case GeomAbs_C1 : return GeomAbs_C0; 
517   case GeomAbs_C2 : return GeomAbs_C1;
518   case GeomAbs_C3 : return GeomAbs_C2;
519   case GeomAbs_CN : return GeomAbs_CN;  
520   default: break;
521   }
522   return GeomAbs_C0; 
523 }
524
525 //=============================================================================
526 //function : IsClosed
527 //purpose  :
528 //=============================================================================
529 Standard_Boolean Bisector_BisecCC::IsClosed() const
530 {
531   if (curve1->IsClosed()) {
532     if (startIntervals.First() == curve1->FirstParameter() &&
533       endIntervals  .Last () == curve1->LastParameter ()    )
534       return Standard_True;
535   }
536   return Standard_False;
537 }
538
539 //=============================================================================
540 //function : IsPeriodic
541 //purpose  :
542 //=============================================================================
543 Standard_Boolean Bisector_BisecCC::IsPeriodic() const
544 {
545   return Standard_False;
546 }
547
548
549 //=============================================================================
550 //function : Curvature
551 //purpose  :
552 //=============================================================================
553 static Standard_Real Curvature (const Handle(Geom2d_Curve)& C,
554   Standard_Real         U,
555   Standard_Real         Tol)
556 {
557   Standard_Real K1; 
558   gp_Vec2d      D1,D2;
559   gp_Pnt2d      P;
560   C->D2(U,P,D1,D2);
561   Standard_Real Norm2 = D1.SquareMagnitude();
562   if (Norm2 < Tol) {
563     K1 = 0.0;
564   }
565   else {
566     K1   = (D1^D2)/(Norm2*sqrt(Norm2));
567   }
568   return K1;
569 }
570
571 //=============================================================================
572 //function : Value
573 //purpose  : CALCULATE THE CURRENT POINT BY ITERATIVE METHOD.
574 //           ----------------------------------------------
575 //           Calculate the current point, the distance from the current point to
576 //           both curves, the parameters on each curve of the projection  
577 //           of the current point.
578 //
579 //method  : - Find start parameter by using <myPolygon>.
580 //          - Calculate parameter U2 on curve C2 solution of H(U,V)= 0
581 //          - P(U) = F(U,U2)
582 //
583 //          or :
584 //                                 ||P2(v0)P1(u)||**2
585 //           F(u,v) = P1(u) - 1/2 *----------------* N(u)
586 //                                (N(u).P2(v0)P1(u))
587 //
588 //           H(u,v) = (Tu.P1(u)P2(v))**2||Tv||**2 - (Tv.P1(u)P2(v))**2||Tu||**2
589 //
590 //=============================================================================
591 gp_Pnt2d Bisector_BisecCC::ValueAndDist (const Standard_Real  U,
592   Standard_Real& U1,
593   Standard_Real& U2,
594   Standard_Real& Dist) const
595 {       
596   gp_Vec2d T;  
597
598   //-----------------------------------------------
599   // is the polygon reduced to a point or empty?
600   //-----------------------------------------------
601   if (myPolygon.Length() <= 1) {
602     return Extension(U,U1,U2,Dist,T);
603   }
604
605   //-----------------------------------------------
606   // test U out of the limits of the polygon.
607   //-----------------------------------------------
608   if (U  <  myPolygon.First().ParamOnBis()) {
609     return Extension(U,U1,U2,Dist,T);
610   }
611   if (U >  myPolygon.Last().ParamOnBis()) {
612     return Extension(U,U1,U2,Dist,T);
613   }
614
615   //-------------------------------------------------------
616   // Find start parameter by using <myPolygon>.
617   //-------------------------------------------------------
618   Standard_Integer IntervalIndex = myPolygon.Interval(U);
619   Standard_Real    UMin = myPolygon.Value(IntervalIndex    ).ParamOnBis();
620   Standard_Real    UMax = myPolygon.Value(IntervalIndex + 1).ParamOnBis();
621   Standard_Real    VMin = myPolygon.Value(IntervalIndex    ).ParamOnC2();
622   Standard_Real    VMax = myPolygon.Value(IntervalIndex + 1).ParamOnC2(); 
623   Standard_Real    Alpha,VInit;
624
625   if (Abs(UMax - UMin) < gp::Resolution()) {
626     VInit = VMin;
627   }
628   else {
629     Alpha = (U - UMin)/(UMax - UMin);
630     VInit = VMin + Alpha*(VMax - VMin);
631   }
632
633   U1 = LinkBisCurve(U);
634   Standard_Real   VTemp = Min(VMin,VMax); 
635   VMax = Max(VMin,VMax); VMin = VTemp;
636   Standard_Boolean Valid = Standard_True; 
637   //---------------------------------------------------------------
638   // Calculate parameter U2 on curve C2 solution of H(u,v)=0     
639   //---------------------------------------------------------------
640   gp_Pnt2d P1;
641   gp_Vec2d T1;
642   Standard_Real    EpsH     = 1.E-9;
643   Standard_Real    EpsH100  = 1.E-7; 
644   curve1->D1 (U1,P1,T1);
645   gp_Vec2d N1(T1.Y(), - T1.X());
646
647   if ((VMax - VMin) < Precision::PConfusion()) {
648     U2 = VInit;
649   }
650   else {
651     Bisector_FunctionH H  (curve2,P1,sign1*sign2*T1);
652     Standard_Real FInit;
653     H.Value(VInit,FInit);
654     if (Abs(FInit) < EpsH) {
655       U2 = VInit;
656     }
657     else
658     {
659       math_BissecNewton aNewSolution(EpsH);
660       aNewSolution.Perform(H, VMin - EpsH100, VMax + EpsH100, 10);
661
662       if (aNewSolution.IsDone())
663       {
664         U2 = aNewSolution.Root();
665       }
666       else
667       {
668         math_FunctionRoot SolRoot (H,VInit,EpsH,VMin - EpsH100,VMax + EpsH100);
669         
670         if (SolRoot.IsDone())
671           U2    = SolRoot.Root();
672         else
673           Valid = Standard_False;
674       }
675     }
676   }
677
678   gp_Pnt2d PBis = pointStart;
679   //----------------
680   // P(U) = F(U1,U2)
681   //----------------
682   if (Valid) {
683     gp_Pnt2d P2 = curve2->Value(U2);
684     gp_Vec2d P2P1(P1.X() - P2.X(),P1.Y() - P2.Y());
685     Standard_Real SquareP2P1 = P2P1.SquareMagnitude();
686     Standard_Real N1P2P1     = N1.Dot(P2P1);
687     const Standard_Real anEps = Epsilon(1);
688
689     if (P1.IsEqual(P2,Precision::Confusion())) {
690       PBis = P1 ;
691       Dist = 0.0;
692     }
693     else if (N1P2P1*sign1 < anEps) {
694       Valid = Standard_False;
695     }
696     else {      
697       PBis = P1.Translated(- (0.5*SquareP2P1/N1P2P1)*N1);
698       Dist = P1.SquareDistance(PBis);
699     }
700   }
701
702   //----------------------------------------------------------------
703   // If the point is not valid 
704   // calculate by intersection.
705   //----------------------------------------------------------------
706   if (!Valid) {
707     //--------------------------------------------------------------------
708     // Construction of the bisectrice point curve and of the straight line passing
709     // by P1 and carried by the normal. curve2 is limited by VMin and VMax.
710     //--------------------------------------------------------------------
711     Standard_Real DMin = Precision::Infinite();
712     gp_Pnt2d      P;
713     Handle(Bisector_BisecPC) BisPC 
714       = new Bisector_BisecPC(curve2, P1, sign2, VMin, VMax);
715     Handle(Geom2d_Line)      NorLi = new Geom2d_Line (P1,N1);
716
717     Geom2dAdaptor_Curve ABisPC(BisPC);
718     Geom2dAdaptor_Curve ANorLi(NorLi);    
719     //-------------------------------------------------------------------------
720     Geom2dInt_GInter  Intersect(ABisPC,ANorLi,
721       Precision::Confusion(),Precision::Confusion());
722     //-------------------------------------------------------------------------
723
724     if (Intersect.IsDone() && !Intersect.IsEmpty()) {
725       for (Standard_Integer i = 1; i <= Intersect.NbPoints(); i++) {
726         if (Intersect.Point(i).ParamOnSecond()*sign1 < Precision::PConfusion()) {
727           P  = Intersect.Point(i).Value();
728           if (P.SquareDistance(P1) < DMin) {
729             DMin = P.SquareDistance(P1);
730             PBis = P;
731             U2   = BisPC->LinkBisCurve(Intersect.Point(i).ParamOnFirst());
732             Dist = DMin;
733           }
734         }
735       }
736     }
737   } 
738   return   PBis;  
739 }
740
741 //=============================================================================
742 //function : ValueByInt
743 //purpose  : CALCULATE THE CURRENT POINT BY INTERSECTION.
744 //           -------------------------------------------
745 //           Calculate the current point, the distance from the current point
746 //           to two curves, the parameters on each curve of the projection of the 
747 //           current point.
748 //           the current point with parameter U is the intersection of the 
749 //           bissectrice point curve (P1,curve2) and of the straight line 
750 //           passing through P1 of director vector N1.
751 //           P1 is the current point of parameter U on curve1 and N1 the  
752 //           normal at this point.
753 //=============================================================================
754 gp_Pnt2d Bisector_BisecCC::ValueByInt (const Standard_Real  U,
755   Standard_Real& U1,
756   Standard_Real& U2,
757   Standard_Real& Dist) const
758 {
759   //------------------------------------------------------------------
760   // Return point, tangent, normal on C1 at parameter U.
761   //-------------------------------------------------------------------
762   U1 = LinkBisCurve(U);
763
764   gp_Pnt2d P1,P2,P,PSol;
765   gp_Vec2d Tan1,Tan2;
766   curve1->D1(U1,P1,Tan1);
767   gp_Vec2d N1(  Tan1.Y(), - Tan1.X());
768
769   //--------------------------------------------------------------------------
770   // test confusion of P1 with extremity of curve2.
771   //--------------------------------------------------------------------------
772   if (P1.Distance(curve2->Value(curve2->FirstParameter())) < Precision::Confusion()) {
773     U2   = curve2->FirstParameter();
774     curve2->D1(U2,P2,Tan2);       
775     if ( isConvex1 && isConvex2 ) {
776       Dist = 0.;
777       return P1;
778     }
779     if (! Tan1.IsParallel(Tan2,Precision::Angular())) {
780       Dist = 0.;
781       return P1;
782     }
783   }
784   if (P1.Distance(curve2->Value(curve2->LastParameter())) < Precision::Confusion()) {
785     U2   = curve2->LastParameter();  
786     curve2->D1(U2,P2,Tan2);    
787     if ( isConvex1 && isConvex2 ) {
788       Dist = 0.;
789       return P1;
790     }
791     if (! Tan1.IsParallel(Tan2,Precision::Angular())) {  
792       Dist = 0.;
793       return P1;
794     }
795   }
796
797   Standard_Boolean YaSol = Standard_False;
798   Standard_Real    DMin  = Precision::Infinite();
799   Standard_Real    USol;
800   Standard_Real    EpsMax = 1.E-6;
801   Standard_Real    EpsX;
802   Standard_Real    EpsH   = 1.E-8;
803   Standard_Real    DistPP1;
804   Standard_Integer NbSamples =20;
805   Standard_Real    UFirstOnC2 = curve2->FirstParameter();
806   Standard_Real    ULastOnC2  = curve2->LastParameter();
807
808   if (!myPolygon.IsEmpty()){
809     if (sign1 == sign2) { ULastOnC2  = myPolygon.Last().ParamOnC2();}
810     else                { UFirstOnC2 = myPolygon.Last().ParamOnC2();} 
811   }
812
813   if (Abs(ULastOnC2 - UFirstOnC2) < Precision::PConfusion()/100.) {
814     Dist = Precision::Infinite();
815     return P1;
816   }
817
818   DiscretPar(Abs(ULastOnC2 - UFirstOnC2),EpsH,EpsMax,2,20,EpsX,NbSamples);
819
820   Bisector_FunctionH H  (curve2,P1,sign1*sign2*Tan1);
821   math_FunctionRoots SolRoot (H,
822     UFirstOnC2,
823     ULastOnC2 ,
824     NbSamples,
825     EpsX,EpsH,EpsH);
826   if (SolRoot.IsDone()) {
827     for (Standard_Integer j = 1; j <= SolRoot.NbSolutions(); j++) {
828       USol        = SolRoot.Value(j);
829       gp_Pnt2d P2Curve2 = curve2->Value(USol);
830       gp_Vec2d P2P1(P1.X() - P2Curve2.X(),P1.Y() - P2Curve2.Y());
831       Standard_Real SquareP2P1 = P2P1.SquareMagnitude();
832       Standard_Real N1P2P1     = N1.Dot(P2P1);
833
834       // Test if the solution is at the proper side of the curves.
835       if (N1P2P1*sign1 > 0 ) {
836         P       = P1.Translated(- (0.5*SquareP2P1/N1P2P1)*N1);
837         DistPP1 = P1.SquareDistance(P);
838         if (DistPP1  < DMin) {
839           DMin  = DistPP1;
840           PSol  = P;
841           U2    = USol;
842           YaSol = Standard_True;
843         }
844       }
845     }
846   }
847
848   /*
849   if (!YaSol) {
850   //--------------------------------------------------------------------
851   // Construction de la bisectrice point courbe et de la droite passant
852   // par P1 et portee par la normale.
853   //--------------------------------------------------------------------
854   Handle(Bisector_BisecPC) BisPC 
855   = new Bisector_BisecPC(curve2,P1,sign2,2*distMax);
856   //-------------------------------
857   // Test si la bissectrice existe.
858   //-------------------------------
859   if (BisPC->IsEmpty()) {
860   Dist = Precision::Infinite();
861   PSol = P1;
862   return PSol;
863   }
864
865   Handle(Geom2d_Line)      NorLi  = new Geom2d_Line (P1,N1);
866   Geom2dAdaptor_Curve      NorLiAd;
867   if (sign1 < 0.) {NorLiAd.Load(NorLi,0.       ,distMax);}
868   else            {NorLiAd.Load(NorLi,- distMax,0.     );}
869
870   //-------------------------------------------------------------------------
871   Geom2dInt_GInter  Intersect(BisPC,NorLiAd,
872   Precision::Confusion(),Precision::Confusion());
873   //-------------------------------------------------------------------------
874   if (Intersect.IsDone() && !Intersect.IsEmpty()) {
875   for (Standard_Integer i = 1; i <= Intersect.NbPoints(); i++) {
876   if (Intersect.Point(i).ParamOnSecond()*sign1< Precision::PConfusion()) {
877   P       = Intersect.Point(i).Value();
878   DistPP1 = P.SquareDistance(P1);
879   if (DistPP1  < DMin) {
880   DMin  = DistPP1;
881   PSol  = P;
882   U2   = Intersect.Point(i).ParamOnFirst();
883   YaSol = Standard_True;
884   }
885   }
886   }
887   }
888   }
889   */
890
891   if (YaSol) {
892     Dist = DMin;
893     //--------------------------------------------------------------
894     // Point found => Test curve distance + Angular Test 
895     //---------------------------------------------------------------
896     P2 = curve2->Value(U2);
897     gp_Vec2d PP1(P1.X() - PSol.X(),P1.Y() - PSol.Y());
898     gp_Vec2d PP2(P2.X() - PSol.X(),P2.Y() - PSol.Y());
899
900     //-----------------------------------------------
901     // Dist = product of norms = distance at the square.
902     //-----------------------------------------------
903     if (PP1.Dot(PP2) > (1. - Precision::Angular())*Dist) {
904       YaSol = Standard_False;
905     }
906     else {
907       if ( !isConvex1 ) {
908         Standard_Real K1 = Curvature(curve1,U1,Precision::Confusion());
909         if (K1 != 0.) {
910           if (Dist > 1/(K1*K1)) YaSol = Standard_False;
911         }
912       }
913       if (YaSol) {
914         if ( !isConvex2 ) {
915           Standard_Real K2 = Curvature(curve2,U2,Precision::Confusion());
916           if (K2 != 0.) {
917             if (Dist > 1/(K2*K2)) YaSol = Standard_False;
918           }
919         }
920       }
921     }
922   }
923   if (!YaSol) { 
924     Dist = Precision::Infinite();
925     PSol = P1;
926   }
927   return PSol;
928 }
929
930 //=============================================================================
931 //function : D0
932 //purpose : 
933 //=============================================================================
934 void Bisector_BisecCC::D0(const Standard_Real     U,
935   gp_Pnt2d&         P) const
936 {
937   Standard_Real U1,U2,Dist;
938
939   P = ValueAndDist(U,U1,U2,Dist);
940 }
941
942 //=============================================================================
943 //function : D1
944 //purpose : 
945 //=============================================================================
946 void Bisector_BisecCC::D1(const Standard_Real     U,
947   gp_Pnt2d&         P,
948   gp_Vec2d&         V ) const
949 {
950   V.SetCoord(0.,0.);
951   gp_Vec2d V2,V3;
952   Values(U,1,P,V,V2,V3);
953
954
955 //=============================================================================
956 //function : D2
957 //purpose : 
958 //=============================================================================
959 void Bisector_BisecCC::D2(const Standard_Real     U,
960   gp_Pnt2d&         P,
961   gp_Vec2d&         V1,
962   gp_Vec2d&         V2) const
963 {
964   V1.SetCoord(0.,0.);
965   V2.SetCoord(0.,0.);
966   gp_Vec2d V3;
967   Values(U,2,P,V1,V2,V3);
968 }
969
970 //=============================================================================
971 //function : D3
972 //purpose : 
973 //=============================================================================
974 void Bisector_BisecCC::D3(const Standard_Real     U,
975   gp_Pnt2d&         P,
976   gp_Vec2d&         V1,
977   gp_Vec2d&         V2,
978   gp_Vec2d&         V3) const
979 {
980   V1.SetCoord(0.,0.);
981   V2.SetCoord(0.,0.);
982   V3.SetCoord(0.,0.);
983   Values(U,3,P,V1,V2,V3);
984 }
985
986 //=============================================================================
987 //function : DN
988 //purpose : 
989 //=============================================================================
990 gp_Vec2d Bisector_BisecCC::DN(const Standard_Real    U,
991   const Standard_Integer N) const
992 {
993   gp_Pnt2d P;
994   gp_Vec2d V1(0.,0.);
995   gp_Vec2d V2(0.,0.);
996   gp_Vec2d V3(0.,0.);
997   Values (U,N,P,V1,V2,V3);
998   switch (N) {
999   case 1 : return V1;
1000   case 2 : return V2;
1001   case 3 : return V3;
1002   default: {
1003     throw Standard_NotImplemented();
1004            }
1005   }
1006 }
1007
1008 //=============================================================================
1009 //function : Values
1010 // purpose : the curve can be described by the following equations:
1011 //
1012 //           B(u) = F(u,v0)
1013 //           where v0 = Phi(u) is given by H (u,v) = 0.
1014 //   
1015 //           with :
1016 //                                 ||P2(v0)P1(u)||**2
1017 //           F(u,v) = P1(u) - 1/2 *----------------* N(u)
1018 //                                (N(u).P2(v0)P1(u))
1019 //
1020 //           H(u,v) = (Tu.P1(u)P2(v))**2||Tv||**2 - (Tv.P1(u)P2(v))**2||Tu||**2
1021 //
1022 //           => dB(u)/du = dF/du + dF/dv(- dH/du:dH/dv)
1023 //
1024 //           Note : tangent to the bisectrice is bissectrice at 
1025 //                      tangents  T1(u) and T2(v0)
1026 //
1027 //=============================================================================
1028 void  Bisector_BisecCC::Values (const Standard_Real     U,
1029   const Standard_Integer  N,
1030   gp_Pnt2d&         P,
1031   gp_Vec2d&         V1,
1032   gp_Vec2d&         V2,
1033   gp_Vec2d&         V3) const
1034 {
1035   V1 = gp_Vec2d(0.,0.);
1036   V2 = gp_Vec2d(0.,0.);
1037   V3 = gp_Vec2d(0.,0.);
1038   //-------------------------------------------------------------------------
1039   // Calculate the current point on the bisectrice and the parameters on each  
1040   // curve.
1041   //-------------------------------------------------------------------------
1042   Standard_Real U0,V0,Dist;  
1043
1044   //-----------------------------------------------
1045   // is the polygon reduced to a point or empty?
1046   //-----------------------------------------------
1047   if (myPolygon.Length() <= 1) {
1048     P = Extension(U,U0,V0,Dist,V1);
1049   }
1050   if (U  <  myPolygon.First().ParamOnBis()) {
1051     P = Extension(U,U0,V0,Dist,V1);
1052     return;
1053   }
1054   if (U >  myPolygon.Last().ParamOnBis()) {
1055     P = Extension(U,U0,V0,Dist,V1);
1056     return;
1057   }
1058   P = ValueAndDist(U,U0,V0,Dist);
1059
1060   if (N == 0) return;
1061   //------------------------------------------------------------------
1062   // Return point, tangent, normal to C1 by parameter U0.
1063   //-------------------------------------------------------------------
1064   gp_Pnt2d P1  ;      // point on C1. 
1065   gp_Vec2d Tu  ;      // tangent to C1 by U0.
1066   gp_Vec2d Tuu ;      // second derivative to C1 by U0.
1067   curve1->D2(U0,P1,Tu,Tuu);  
1068   gp_Vec2d Nor( - Tu .Y() , Tu .X());           // Normal by U0.
1069   gp_Vec2d Nu ( - Tuu.Y() , Tuu.X());           // derivative of the normal by U0.
1070
1071   //-------------------------------------------------------------------
1072   // Return point, tangent, normale to C2 by parameter V0.
1073   //-------------------------------------------------------------------
1074   gp_Pnt2d P2  ;      // point on C2. 
1075   gp_Vec2d Tv  ;      // tangent to C2 by V.
1076   gp_Vec2d Tvv ;      // second derivative to C2 by V.
1077   curve2->D2(V0,P2,Tv,Tvv);  
1078
1079   gp_Vec2d PuPv(P2.X() - P1.X(), P2.Y() - P1.Y());
1080
1081   //-----------------------------
1082   // Calculate dH/du and dH/dv.
1083   //-----------------------------
1084   Standard_Real TuTu,TvTv,TuTv;
1085   Standard_Real TuPuPv,TvPuPv  ;
1086   Standard_Real TuuPuPv,TuTuu  ;
1087   Standard_Real TvvPuPv,TvTvv  ;
1088
1089   TuTu    = Tu.Dot(Tu)    ; TvTv   = Tv.Dot(Tv)  ; TuTv = Tu.Dot(Tv);
1090   TuPuPv  = Tu.Dot(PuPv)  ; TvPuPv = Tv.Dot(PuPv);
1091   TuuPuPv = Tuu.Dot(PuPv) ; TuTuu  = Tu.Dot(Tuu) ;
1092   TvvPuPv = Tvv.Dot(PuPv) ; TvTvv  = Tv.Dot(Tvv) ;
1093
1094   Standard_Real dHdu = 2*(TuPuPv*(TuuPuPv - TuTu)*TvTv +
1095     TvPuPv*TuTv*TuTu  -TuTuu*TvPuPv*TvPuPv);
1096   Standard_Real dHdv = 2*(TuPuPv*TuTv*TvTv + TvTvv*TuPuPv*TuPuPv -
1097     TvPuPv*(TvvPuPv + TvTv)*TuTu);
1098
1099   //-----------------------------
1100   // Calculate dF/du and dF/dv.
1101   //-----------------------------
1102   Standard_Real NorPuPv,NuPuPv,NorTv;
1103   Standard_Real A,B,dAdu,dAdv,dBdu,dBdv,BB;
1104
1105   NorPuPv = Nor.Dot(PuPv); 
1106   NuPuPv  = Nu .Dot(PuPv);  
1107   NorTv   = Nor.Dot(Tv)  ;
1108
1109   A    = 0.5*PuPv.SquareMagnitude();
1110   B    = - NorPuPv;
1111   BB   = B*B;
1112   dAdu = - TuPuPv;
1113   dBdu = - NuPuPv ;
1114   dAdv = TvPuPv;
1115   dBdv = - NorTv;
1116
1117   //---------------------------------------
1118   // F(u,v) = Pu - (A(u,v)/B(u,v))*Nor(u)
1119   //----------------------------------------
1120   if (BB < gp::Resolution()) {
1121     V1 = Tu.Normalized() + Tv.Normalized();
1122     V1 = 0.5*Tu.SquareMagnitude()*V1;
1123   }
1124   else {
1125     gp_Vec2d dFdu = Tu - (dAdu/B - dBdu*A/BB)*Nor - (A/B)*Nu; 
1126     gp_Vec2d dFdv = ( - dAdv/B + dBdv*A/BB)*Nor ;
1127
1128     if (Abs(dHdv) > gp::Resolution()) {
1129       V1 = dFdu + dFdv*( - dHdu / dHdv );
1130     }
1131     else {
1132       V1 = Tu;
1133     }
1134   }
1135   if (N == 1) return;
1136 }
1137
1138 //=============================================================================
1139 //function : Extension
1140 // purpose : Calculate the current point on the  extensions 
1141 //           by tangence of the curve.
1142 //============================================================================
1143 gp_Pnt2d Bisector_BisecCC::Extension (const Standard_Real  U,
1144   Standard_Real& U1,
1145   Standard_Real& U2,
1146   Standard_Real& Dist,
1147   gp_Vec2d&      T   ) const
1148 {
1149   Bisector_PointOnBis PRef;
1150   gp_Pnt2d            P,P1,P2,PBis;
1151   gp_Vec2d            T1,Tang;
1152   Standard_Real       dU = 0.;
1153   Standard_Boolean    ExtensionTangent = Standard_False;
1154
1155   if (myPolygon.Length() == 0) {
1156     //---------------------------------------------
1157     // Empty Polygon => segment (pointStart,pointEnd)
1158     //---------------------------------------------
1159     dU = U - startIntervals.First();
1160     P  = pointStart;
1161     P1 = pointEnd;
1162     U1 = curve1->LastParameter();
1163     if (sign1 == sign2) { U2 = curve2->FirstParameter();}
1164     else                { U2 = curve2->LastParameter() ;}
1165     Tang.SetCoord(P1.X() - P.X(),P1.Y() - P.Y());
1166   }
1167   else if (U  < myPolygon.First().ParamOnBis()) {
1168     PRef = myPolygon.First();    
1169     P    = PRef.Point();    
1170     dU   = U - PRef.ParamOnBis();
1171     if (extensionStart) {
1172       //------------------------------------------------------------
1173       // extension = segment (pointstart, first point of the polygon.)
1174       //------------------------------------------------------------
1175       P1 = pointStart;
1176       U1 = curve1->FirstParameter();
1177       if (sign1 == sign2) { U2 = curve2->LastParameter();}
1178       else                { U2 = curve2->FirstParameter();} 
1179       Tang.SetCoord(P.X() - P1.X(),P.Y() - P1.Y());
1180     }
1181     else {
1182       ExtensionTangent = Standard_True;
1183     }
1184   }
1185   else if (U >  myPolygon.Last().ParamOnBis()) {
1186     PRef = myPolygon.Last();    
1187     P    = PRef.Point();    
1188     dU   = U - PRef.ParamOnBis();
1189     if (extensionEnd) {      
1190       //------------------------------------------------------------
1191       // extension = segment (last point of the polygon.pointEnd)
1192       //------------------------------------------------------------
1193       P1 = pointEnd;
1194       U1 = curve1->LastParameter();
1195       if (sign1 == sign2) { U2 = curve2->LastParameter();}
1196       else                { U2 = curve2->FirstParameter();} 
1197       Tang.SetCoord(P1.X() - P.X(),P1.Y() - P.Y());
1198     }
1199     else {
1200       ExtensionTangent = Standard_True;     
1201     }
1202   }
1203
1204   if (ExtensionTangent) {
1205     //-----------------------------------------------------------
1206     // If the la curve has no a extension, it is extended by tangency      
1207     //------------------------------------------------------------
1208     U1 = PRef.ParamOnC1();
1209     U2 = PRef.ParamOnC2(); 
1210     P2 = curve2->Value(U2);
1211     curve1->D1(U1,P1,T1);
1212     Tang.SetCoord(2*P.X() - P1.X() - P2.X(), 2*P.Y() - P1.Y() - P2.Y());
1213     if (Tang.Magnitude() < Precision::Confusion()) {
1214       Tang = T1;
1215     }
1216     if (T1.Dot(Tang) < 0.) Tang = - Tang;
1217   } 
1218
1219   T = Tang.Normalized();
1220   PBis.SetCoord(P.X() + dU*T.X(),P.Y() + dU*T.Y());        
1221   Dist = P1.Distance(PBis);
1222   return PBis;
1223 }
1224
1225 //=============================================================================
1226 //function : PointByInt
1227 // purpose : 
1228 //=============================================================================
1229 static Standard_Boolean PointByInt(const Handle(Geom2d_Curve)& CA,
1230   const Handle(Geom2d_Curve)& CB,
1231   const Standard_Real         SignA,
1232   const Standard_Real         SignB,
1233   const Standard_Real         UOnA,
1234   Standard_Real&        UOnB,
1235   Standard_Real&        Dist)
1236 {  
1237   //------------------------------------------------------------------
1238   // Return point,tangent, normal on CA with parameter UOnA.
1239   //-------------------------------------------------------------------
1240   gp_Pnt2d P1,P2,P,PSol;
1241   gp_Vec2d Tan1,Tan2;
1242   Standard_Boolean IsConvexA = Bisector::IsConvex(CA,SignA);
1243   Standard_Boolean IsConvexB = Bisector::IsConvex(CB,SignB);
1244
1245   CA->D1(UOnA,P1,Tan1);
1246   gp_Vec2d N1(Tan1.Y(), - Tan1.X());
1247
1248   //--------------------------------------------------------------------------
1249   // test of confusion of P1 with extremity of curve2.
1250   //--------------------------------------------------------------------------
1251   if (P1.Distance(CB->Value(CB->FirstParameter())) < Precision::Confusion()) {
1252     UOnB = CB->FirstParameter();
1253     CB->D1(UOnB,P2,Tan2);         
1254     if ( IsConvexA && IsConvexB ) {
1255       Dist = 0.;
1256       return Standard_True;
1257     }
1258     if (! Tan1.IsParallel(Tan2,Precision::Angular())) {
1259       Dist = 0.;
1260       return Standard_False;
1261     }
1262   }
1263   if (P1.Distance(CB->Value(CB->LastParameter())) < Precision::Confusion()) {
1264     UOnB   = CB->LastParameter();  
1265     CB->D1(UOnB,P2,Tan2);    
1266     if ( IsConvexA && IsConvexB ) {
1267       Dist = 0.;
1268       return Standard_True;
1269     }
1270     if (! Tan1.IsParallel(Tan2,Precision::Angular())) {  
1271       Dist = 0.;
1272       return Standard_False;
1273     }
1274   }
1275
1276   Standard_Real    DMin = Precision::Infinite();
1277   Standard_Real    UPC;
1278   Standard_Boolean YaSol = Standard_False; 
1279   //--------------------------------------------------------------------
1280   // Construction of the bisectrice point curve and of the straight line passing
1281   // through P1 and carried by the normal. 
1282   //--------------------------------------------------------------------
1283   Handle(Bisector_BisecPC) BisPC 
1284     = new Bisector_BisecPC(CB,P1,SignB );
1285   //-------------------------------
1286   // Test if the bissectrice exists.
1287   //-------------------------------  
1288   if (BisPC->IsEmpty()) {
1289     Dist = Precision::Infinite();
1290     PSol = P1;
1291     return Standard_False;
1292   }
1293
1294   Handle(Geom2d_Line)      NorLi = new Geom2d_Line (P1,N1);
1295
1296   Geom2dAdaptor_Curve ABisPC(BisPC);
1297   Geom2dAdaptor_Curve ANorLi(NorLi);    
1298   //-------------------------------------------------------------------------
1299   Geom2dInt_GInter  Intersect(ABisPC,ANorLi,
1300     Precision::Confusion(),Precision::Confusion());
1301   //-------------------------------------------------------------------------
1302
1303   if (Intersect.IsDone() && !Intersect.IsEmpty()) {
1304     for (Standard_Integer i = 1; i <= Intersect.NbPoints(); i++) {
1305       if (Intersect.Point(i).ParamOnSecond()*SignA < Precision::PConfusion()) {
1306         P  = Intersect.Point(i).Value();
1307         if (P.SquareDistance(P1) < DMin) {
1308           DMin  = P.SquareDistance(P1);
1309           PSol  = P;
1310           UPC   = Intersect.Point(i).ParamOnFirst();
1311           UOnB  = BisPC->LinkBisCurve(UPC);
1312           Dist  = DMin;
1313           YaSol = Standard_True;
1314         }
1315       }
1316     }
1317   }  
1318   if (YaSol) {    
1319     //--------------------------------------------------------------
1320     // Point found => Test distance curvature + Angular test 
1321     //---------------------------------------------------------------
1322     P2 = CB->Value(UOnB);
1323     if(P1.SquareDistance(PSol) < 1.e-32)
1324     {
1325       YaSol = Standard_False;
1326       return YaSol;
1327     }
1328     if(P2.SquareDistance(PSol) < 1.e-32)
1329     {
1330       YaSol = Standard_False;
1331       return YaSol;
1332     }
1333
1334     gp_Dir2d PP1Unit(P1.X() - PSol.X(),P1.Y() - PSol.Y());
1335     gp_Dir2d PP2Unit(P2.X() - PSol.X(),P2.Y() - PSol.Y());
1336
1337     if (PP1Unit*PP2Unit > 1. - Precision::Angular()) {
1338       YaSol = Standard_False;
1339     }
1340     else {
1341       Dist = sqrt(Dist);
1342       if ( !IsConvexA ) {
1343         Standard_Real K1 = Curvature(CA,UOnA,Precision::Confusion());
1344         if (K1 != 0.) {
1345           if (Dist > Abs(1/K1)) YaSol = Standard_False;
1346         }
1347       }
1348       if (YaSol) {
1349         if ( !IsConvexB ) {
1350           Standard_Real K2 = Curvature(CB,UOnB,Precision::Confusion());
1351           if (K2 != 0.) {
1352             if (Dist > Abs(1/K2)) YaSol = Standard_False;
1353           }
1354         }
1355       }
1356     }
1357   }
1358   return YaSol;
1359 }
1360
1361 //=============================================================================
1362 //function : SupLastParameter
1363 // purpose : 
1364 //=============================================================================
1365 void Bisector_BisecCC::SupLastParameter()
1366 {
1367   endIntervals.Append(curve1->LastParameter());
1368   // ----------------------------------------------------------------------
1369   // Calculate parameter on curve1 associated to one or the other of the extremities
1370   // of curve2 following the values of sign1 and sign2.
1371   // the bissectrice is limited by the obtained parameters.
1372   //------------------------------------------------------------------------
1373   Standard_Real    UOnC1,UOnC2,Dist;
1374   if (sign1 == sign2) {
1375     UOnC2 =  curve2->FirstParameter();
1376   }
1377   else {
1378     UOnC2 = curve2->LastParameter();
1379   }
1380   Standard_Boolean YaSol  = PointByInt(curve2,curve1,sign2,sign1,UOnC2,UOnC1,Dist);
1381   if (YaSol) {
1382     if (UOnC1 > startIntervals.First() && UOnC1 < endIntervals.Last()) {
1383       endIntervals.SetValue(1,UOnC1);
1384     }
1385   }
1386 }
1387
1388 //=============================================================================
1389 //function : Curve
1390 // purpose : 
1391 //=============================================================================
1392 Handle(Geom2d_Curve) Bisector_BisecCC::Curve(const Standard_Integer I) const 
1393 {
1394   if      (I == 1) return curve1;
1395   else if (I == 2) return curve2;
1396   else             throw Standard_OutOfRange();
1397 }
1398
1399 //=============================================================================
1400 //function : LinkBisCurve
1401 //purpose : 
1402 //=============================================================================
1403 Standard_Real Bisector_BisecCC::LinkBisCurve(const Standard_Real U) const 
1404 {
1405   return (U - shiftParameter);
1406
1407
1408 //=============================================================================
1409 //function : LinkCurveBis
1410 //purpose : 
1411 //=============================================================================
1412 Standard_Real Bisector_BisecCC::LinkCurveBis(const Standard_Real U) const 
1413 {
1414   return (U + shiftParameter);
1415 }
1416
1417 //=============================================================================
1418 //function : Indent
1419 //purpose  : 
1420 //=============================================================================
1421 static void Indent(const Standard_Integer Offset) {
1422   if (Offset > 0) {
1423     for (Standard_Integer i = 0; i < Offset; i++) {std::cout << " ";}
1424   }
1425 }
1426
1427 //=============================================================================
1428 //function : Polygon
1429 // purpose : 
1430 //=============================================================================
1431 const Bisector_PolyBis&  Bisector_BisecCC::Polygon() const
1432 {
1433   return myPolygon;
1434 }
1435
1436 //==========================================================================
1437 //function : Parameter
1438 //purpose  :
1439 //==========================================================================
1440 Standard_Real Bisector_BisecCC::Parameter(const gp_Pnt2d& P) const
1441 {
1442   Standard_Real UOnCurve;
1443
1444   if (P.IsEqual(Value(FirstParameter()),Precision::Confusion())) {
1445     UOnCurve = FirstParameter();
1446   }
1447   else if (P.IsEqual(Value(LastParameter()),Precision::Confusion())) {
1448     UOnCurve = LastParameter();
1449   }
1450   else
1451   {
1452     ProjOnCurve(P, curve1, UOnCurve);
1453   }
1454
1455   return UOnCurve;
1456 }
1457
1458
1459 //=============================================================================
1460 //function : Dump
1461 // purpose : 
1462 //=============================================================================
1463 //void Bisector_BisecCC::Dump(const Standard_Integer Deep, 
1464 void Bisector_BisecCC::Dump(const Standard_Integer , 
1465   const Standard_Integer Offset) const 
1466 {
1467   Indent (Offset);
1468   std::cout <<"Bisector_BisecCC :"<<std::endl;
1469   Indent (Offset);
1470   //  std::cout <<"Curve1 :"<<curve1<<std::endl;
1471   //  std::cout <<"Curve2 :"<<curve2<<std::endl;
1472   std::cout <<"Sign1  :"<<sign1<<std::endl;
1473   std::cout <<"Sign2  :"<<sign2<<std::endl;
1474
1475   std::cout <<"Number Of Intervals :"<<startIntervals.Length()<<std::endl;
1476   for (Standard_Integer i = 1; i <= startIntervals.Length(); i++) {
1477     std::cout <<"Interval number :"<<i<<"Start :"<<startIntervals.Value(i)
1478       <<"  end :"<<  endIntervals.Value(i)<<std::endl ;
1479   }
1480   std::cout <<"Index Current Interval :"<<currentInterval<<std::endl;
1481 }
1482
1483 //=============================================================================
1484 //function : Curve
1485 // purpose : 
1486 //=============================================================================
1487 void Bisector_BisecCC::Curve(const Standard_Integer      I,
1488   const Handle(Geom2d_Curve)& C)  
1489 {
1490   if      (I == 1) curve1 = C;
1491   else if (I == 2) curve2 = C;
1492   else             throw Standard_OutOfRange();
1493 }
1494
1495 //=============================================================================
1496 //function : Sign
1497 // purpose : 
1498 //=============================================================================
1499 void Bisector_BisecCC::Sign(const Standard_Integer      I,
1500   const Standard_Real         S)
1501 {
1502   if      (I == 1) sign1 = S;
1503   else if (I == 2) sign2 = S;
1504   else             throw Standard_OutOfRange();
1505 }
1506
1507 //=============================================================================
1508 //function : Polygon
1509 // purpose : 
1510 //=============================================================================
1511 void Bisector_BisecCC::Polygon(const Bisector_PolyBis& P)  
1512 {
1513   myPolygon = P;
1514 }
1515
1516 //=============================================================================
1517 //function : DistMax
1518 // purpose : 
1519 //=============================================================================
1520 void Bisector_BisecCC::DistMax(const Standard_Real D)  
1521 {
1522   distMax = D;
1523 }
1524
1525 //=============================================================================
1526 //function : IsConvex
1527 // purpose : 
1528 //=============================================================================
1529 void Bisector_BisecCC::IsConvex(const Standard_Integer     I,
1530   const Standard_Boolean     IsConvex)  
1531 {
1532   if      (I == 1) isConvex1 = IsConvex;
1533   else if (I == 2) isConvex2 = IsConvex;
1534   else             throw Standard_OutOfRange();
1535 }
1536
1537 //=============================================================================
1538 //function : IsEmpty
1539 // purpose : 
1540 //=============================================================================
1541 void Bisector_BisecCC::IsEmpty ( const Standard_Boolean     IsEmpty)  
1542 {
1543   isEmpty = IsEmpty;
1544 }
1545
1546 //=============================================================================
1547 //function : ExtensionStart
1548 // purpose : 
1549 //=============================================================================
1550 void Bisector_BisecCC::ExtensionStart( const Standard_Boolean  ExtensionStart)  
1551 {
1552   extensionStart = ExtensionStart;
1553 }
1554
1555 //=============================================================================
1556 //function : ExtensionEnd
1557 // purpose : 
1558 //=============================================================================
1559 void Bisector_BisecCC::ExtensionEnd( const Standard_Boolean  ExtensionEnd)  
1560 {
1561   extensionEnd = ExtensionEnd;
1562 }
1563
1564 //=============================================================================
1565 //function : PointStart
1566 // purpose : 
1567 //=============================================================================
1568 void Bisector_BisecCC::PointStart( const gp_Pnt2d& Point)  
1569 {
1570   pointStart = Point;
1571 }
1572
1573 //=============================================================================
1574 //function : PointEnd
1575 // purpose : 
1576 //=============================================================================
1577 void Bisector_BisecCC::PointEnd( const gp_Pnt2d& Point)  
1578 {
1579   pointEnd = Point;
1580 }
1581
1582 //=============================================================================
1583 //function : StartIntervals
1584 // purpose : 
1585 //=============================================================================
1586 void Bisector_BisecCC::StartIntervals 
1587   (const TColStd_SequenceOfReal& StartIntervals)  
1588 {
1589   startIntervals = StartIntervals;
1590 }
1591
1592 //=============================================================================
1593 //function : EndIntervals
1594 // purpose : 
1595 //=============================================================================
1596 void Bisector_BisecCC::EndIntervals 
1597   (const TColStd_SequenceOfReal& EndIntervals)  
1598 {
1599   endIntervals = EndIntervals;
1600 }
1601
1602 //=============================================================================
1603 //function : FirstParameter
1604 // purpose : 
1605 //=============================================================================
1606 void Bisector_BisecCC::FirstParameter (const Standard_Real U)  
1607 {
1608   startIntervals.Append(U);
1609 }
1610
1611 //=============================================================================
1612 //function : LastParameter
1613 // purpose : 
1614 //=============================================================================
1615 void Bisector_BisecCC::LastParameter (const Standard_Real U)  
1616 {
1617   endIntervals.Append(U);
1618 }
1619
1620 //=============================================================================
1621 //function : SearchBound
1622 // purpose : 
1623 //=============================================================================
1624 Standard_Real Bisector_BisecCC::SearchBound (const Standard_Real U1,
1625   const Standard_Real U2) const
1626 {
1627   Standard_Real UMid,Dist1,Dist2,DistMid,U11,U22;
1628   Standard_Real UC1,UC2;
1629   gp_Pnt2d PBis,PBisPrec;
1630   Standard_Real TolPnt   = Precision::Confusion();
1631   Standard_Real TolPar   = Precision::PConfusion();
1632   U11 = U1; U22 = U2;
1633   PBisPrec = ValueByInt(U11,UC1,UC2,Dist1);
1634   PBis     = ValueByInt(U22,UC1,UC2,Dist2);
1635
1636   while ((U22 - U11) > TolPar || 
1637     ((Dist1 < Precision::Infinite() && 
1638     Dist2 < Precision::Infinite() &&
1639     !PBis.IsEqual(PBisPrec,TolPnt)))) { 
1640       PBisPrec = PBis;
1641       UMid     = 0.5*( U22 + U11);
1642       PBis     = ValueByInt(UMid,UC1,UC2,DistMid);
1643       if ((Dist1 < Precision::Infinite()) == (DistMid < Precision::Infinite())) {
1644         U11    = UMid;
1645         Dist1 = DistMid;
1646       }
1647       else {
1648         U22    = UMid;
1649         Dist2 = DistMid;
1650       }
1651   }
1652   PBis = ValueByInt(U11,UC1,UC2,Dist1);
1653   if (Dist1 < Precision::Infinite()) {
1654     UMid = U11;
1655   }
1656   else {
1657     UMid = U22;
1658   }
1659   return UMid;
1660 }
1661
1662 //=============================================================================
1663 //function : ProjOnCurve
1664 // purpose :
1665 //=============================================================================
1666 static Standard_Boolean ProjOnCurve (const gp_Pnt2d& P,
1667   const Handle(Geom2d_Curve)& C,
1668   Standard_Real& theParam)
1669 {
1670   //Standard_Real UOnCurve =0.;
1671   theParam = 0.0;
1672   gp_Pnt2d      PF,PL;
1673   gp_Vec2d      TF,TL;
1674
1675   C->D1(C->FirstParameter(),PF,TF);
1676   C->D1(C->LastParameter() ,PL,TL);
1677
1678   if (P.IsEqual(PF ,Precision::Confusion()))
1679   {
1680     theParam =  C->FirstParameter();
1681     return Standard_True;
1682   }
1683
1684   if (P.IsEqual(PL ,Precision::Confusion()))
1685   {
1686     theParam = C->LastParameter();
1687     return Standard_True;
1688   }
1689
1690   gp_Vec2d PPF(PF.X() - P.X(), PF.Y() - P.Y());
1691   TF.Normalize();
1692
1693   if ( Abs (PPF.Dot(TF)) < Precision::Confusion())
1694   {
1695     theParam = C->FirstParameter();
1696     return Standard_True;
1697   }
1698   gp_Vec2d PPL (PL.X() - P.X(), PL.Y() - P.Y());
1699   TL.Normalize();
1700   if ( Abs (PPL.Dot(TL)) < Precision::Confusion())
1701   {
1702     theParam = C->LastParameter();
1703     return Standard_True;
1704   }
1705   Geom2dAPI_ProjectPointOnCurve Proj(P,C,
1706     C->FirstParameter(),
1707     C->LastParameter());
1708   if (Proj.NbPoints() > 0) {
1709     theParam = Proj.LowerDistanceParameter();
1710   }
1711   else {
1712     return Standard_False;
1713   }
1714
1715   return Standard_True;
1716 }
1717
1718 //=============================================================================
1719 //function : TestExtension
1720 // purpose : 
1721 //=============================================================================
1722 static Standard_Boolean TestExtension (const Handle(Geom2d_Curve)& C1,
1723   const Handle(Geom2d_Curve)& C2,
1724   const Standard_Integer      Start_End)
1725 {
1726   gp_Pnt2d         P1,P2;
1727   gp_Vec2d         T1,T2;
1728   Standard_Boolean Test = Standard_False;
1729   if (Start_End == 1) {
1730     C1->D1(C1->FirstParameter(),P1,T1);
1731   }
1732   else {
1733     C1->D1(C1->LastParameter(),P1,T1); 
1734   }
1735   C2->D1(C2->FirstParameter(),P2,T2);
1736   if (P1.IsEqual(P2,Precision::Confusion())) {
1737     T1.Normalize(); T2.Normalize();
1738     if (T1.Dot(T2) > 1.0 - Precision::Confusion()) {
1739       Test = Standard_True;
1740     }
1741   }
1742   else {
1743     C2->D1(C2->LastParameter(),P2,T2);
1744     if (P1.IsEqual(P2,Precision::Confusion())) {
1745       T2.Normalize();
1746       if (T1.Dot(T2) > 1.0 - Precision::Confusion()) {
1747         Test = Standard_True;
1748       }
1749     }
1750   }
1751   return Test;
1752 }
1753
1754 //=============================================================================
1755 //function : ComputePointEnd
1756 // purpose : 
1757 //=============================================================================
1758 void  Bisector_BisecCC::ComputePointEnd  ()
1759 {  
1760   Standard_Real U1,U2;  
1761   Standard_Real KC,RC;
1762   U1 = curve1->FirstParameter();
1763   if (sign1 == sign2) {
1764     U2 = curve2->LastParameter();
1765   }
1766   else {
1767     U2 = curve2->FirstParameter();
1768   }
1769   Standard_Real K1 = Curvature(curve1,U1,Precision::Confusion());
1770   Standard_Real K2 = Curvature(curve2,U2,Precision::Confusion());    
1771   if (!isConvex1 && !isConvex2) {
1772     if (K1 < K2) {KC = K1;} else {KC = K2;}
1773   }
1774   else if (!isConvex1) {KC = K1;}
1775   else                 {KC = K2;}
1776
1777   gp_Pnt2d      PF; 
1778   gp_Vec2d      TF;
1779   curve1->D1(U1,PF,TF);
1780   TF.Normalize();
1781   if (KC != 0.) { RC = Abs(1/KC);}
1782   else          { RC = Precision::Infinite();}
1783   pointEnd.SetCoord(PF.X() - sign1*RC*TF.Y(), PF.Y() + sign1*RC*TF.X());
1784
1785 }
1786
1787 //=============================================================================
1788 //function : DiscretPar
1789 // purpose :
1790 //=============================================================================
1791 static Standard_Boolean DiscretPar(const Standard_Real     DU, 
1792   const Standard_Real     EpsMin, 
1793   const Standard_Real     EpsMax,
1794   const Standard_Integer  NbMin,
1795   const Standard_Integer  NbMax,
1796   Standard_Real&    Eps,
1797   Standard_Integer& Nb) 
1798 {
1799   if (DU <= NbMin*EpsMin) {
1800     Eps = DU/(NbMin + 1) ;
1801     Nb  = NbMin;  
1802     return Standard_False;
1803   }
1804
1805   Eps = Min (EpsMax,DU/NbMax);
1806
1807   if (Eps < EpsMin) {
1808     Eps = EpsMin;
1809     Nb  = Standard_Integer(DU/EpsMin);
1810   }  
1811   else { Nb = NbMax;}
1812
1813   return Standard_True;
1814 }
1815
1816