src/Bisector/Bisector_BisecCC.cxx

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