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