src/Bisector/Bisector_BisecCC.cxx

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