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