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