src/Bisector/Bisector_BisecCC.cxx

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