src/Extrema/Extrema_GExtPC.gxx

   1 // Created on: 1992-10-19
   2 // Created by: Laurent PAINNOT
   3 // Copyright (c) 1992-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 // Modified by cma, Fri Dec 10 18:11:56 1993
  18
  19 #include Extrema_EPC_hxx
  20 #include ThePOnC_hxx
  21 #include TheSequenceOfPOnC_hxx
  22 #include ThePoint_hxx
  23 #include TheVector_hxx
  24 #include TheExtPElC_hxx
  25 #include <StdFail_NotDone.hxx>
  26 #include <Standard_Failure.hxx>
  27 #include <GeomAbs_CurveType.hxx>
  28 #include <Precision.hxx>
  29 #include <ElCLib.hxx>
  30 #include <TColStd_Array1OfReal.hxx>
  31 #include <NCollection_Array1.hxx>
  32
  33
  34 //=======================================================================
  35 //function : Perform
  36 //purpose  :
  37 //=======================================================================
  38 void Extrema_GExtPC::Perform(const ThePoint& P)
  39 {
  40   mySqDist.Clear();
  41   mypoint.Clear();
  42   myismin.Clear();
  43   Standard_Integer i, NbExt, n;
  44   Standard_Real U;
  45   mysample = 17;
  46   Standard_Real t3d = Precision::Confusion();
  47   if (Precision::IsInfinite(myuinf)) mydist1 = RealLast();
  48   else {
  49     Pf = TheCurveTool::Value(*((TheCurve*)myC), myuinf);
  50     mydist1 = P.SquareDistance(Pf);
  51   }
  52
  53   if (Precision::IsInfinite(myusup)) mydist2 = RealLast();
  54   else {
  55     Pl = TheCurveTool::Value(*((TheCurve*)myC), myusup);
  56     mydist2 = P.SquareDistance(Pl);
  57   }
  58
  59   TheCurve & aCurve = *((TheCurve*)myC);
  60
  61   switch(type) {
  62   case GeomAbs_Circle:
  63     {
  64       myExtPElC.Perform(P, TheCurveTool::Circle(aCurve), t3d, myuinf, myusup);
  65     }
  66     break;
  67   case GeomAbs_Ellipse:
  68     {
  69       myExtPElC.Perform(P, TheCurveTool::Ellipse(aCurve), t3d, myuinf, myusup);
  70     }
  71     break;
  72   case GeomAbs_Parabola:
  73     {
  74       myExtPElC.Perform(P, TheCurveTool::Parabola(aCurve), t3d,myuinf,myusup);
  75     }
  76     break;
  77   case GeomAbs_Hyperbola:
  78     {
  79       myExtPElC.Perform(P,TheCurveTool::Hyperbola(aCurve),t3d, myuinf, myusup);
  80     }
  81     break;
  82   case GeomAbs_Line:
  83     {
  84       myExtPElC.Perform(P, TheCurveTool::Line(aCurve), t3d, myuinf, myusup);
  85     }
  86     break;
  87   case GeomAbs_BezierCurve:
  88     {
  89       myintuinf = myuinf;
  90       myintusup = myusup;
  91       mysample = (TheCurveTool::Bezier(aCurve))->NbPoles() * 2;
  92       myExtPC.Initialize(aCurve);
  93       IntervalPerform(P);
  94       return;
  95     }
  96   case GeomAbs_BSplineCurve:
  97     {
  98       const Standard_Integer
  99         aFirstIdx = TheCurveTool::BSpline(aCurve)->FirstUKnotIndex(),
 100         aLastIdx  = TheCurveTool::BSpline(aCurve)->LastUKnotIndex();
 101       // const reference can not be used due to implementation of BRep_Adaptor.
 102       TColStd_Array1OfReal aKnots(aFirstIdx, aLastIdx);
 103       TheCurveTool::BSpline(aCurve)->Knots(aKnots);
 104
 105       // Workaround to work with:
 106       // blend, where knots may be moved from param space.
 107       Standard_Real aPeriodJump = 0.0;
 108       if (TheCurveTool::IsPeriodic(aCurve))
 109       {
 110         Standard_Integer aPeriodShift =
 111           Standard_Integer ((myuinf - aKnots(aFirstIdx)) / TheCurveTool::Period(aCurve));
 112
 113         if (myuinf < aKnots(aFirstIdx))
 114           aPeriodShift--;
 115
 116         aPeriodJump = TheCurveTool::Period(aCurve) * aPeriodShift;
 117       }
 118
 119       Standard_Integer anIdx;
 120
 121       // Find first and last used knot
 122       Standard_Integer aFirstUsedKnot = aFirstIdx,
 123                         aLastUsedKnot = aLastIdx;
 124       for(anIdx = aFirstIdx; anIdx <= aLastIdx; anIdx++)
 125       {
 126         Standard_Real aKnot = aKnots(anIdx) + aPeriodJump;
 127         if (myuinf >= aKnot)
 128           aFirstUsedKnot = anIdx;
 129         else
 130           break;
 131
 132       }
 133       for(anIdx = aLastIdx; anIdx >= aFirstIdx; anIdx--)
 134       {
 135         Standard_Real aKnot = aKnots(anIdx) + aPeriodJump;
 136         if (myusup <= aKnot)
 137           aLastUsedKnot = anIdx;
 138         else
 139           break;
 140       }
 141
 142       mysample = (TheCurveTool::BSpline(aCurve))->Degree() + 1;
 143
 144       // Fill sample points.
 145       Standard_Integer aValIdx = 1;
 146       NCollection_Array1<Standard_Real> aVal(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
 147       NCollection_Array1<Standard_Real> aParam(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
 148       for(anIdx = aFirstUsedKnot; anIdx < aLastUsedKnot; anIdx++)
 149       {
 150         Standard_Real aF = aKnots(anIdx) + aPeriodJump,
 151                       aL = aKnots(anIdx + 1) + aPeriodJump;
 152
 153         if (anIdx == aFirstUsedKnot)
 154           aF = myuinf;
 155         if (anIdx == aLastUsedKnot - 1)
 156           aL = myusup;
 157
 158         Standard_Real aStep = (aL - aF) / mysample;
 159         for(Standard_Integer aPntIdx = 0; aPntIdx < mysample; aPntIdx++)
 160         {
 161           Standard_Real aCurrentParam = aF + aStep * aPntIdx;
 162           aVal(aValIdx) = TheCurveTool::Value(aCurve, aCurrentParam).SquareDistance(P);
 163           aParam(aValIdx) = aCurrentParam;
 164           aValIdx++;
 165         }
 166       }
 167       // Fill last point.
 168       aVal(aValIdx) = TheCurveTool::Value(aCurve, myusup).SquareDistance(P);
 169       aParam(aValIdx) = myusup;
 170
 171        myExtPC.Initialize(aCurve);
 172
 173       // Find extremas.
 174       for(anIdx = aVal.Lower() + 1; anIdx < aVal.Upper(); anIdx++)
 175       {
 176         if (aVal(anIdx) <= Precision::SquareConfusion())
 177         {
 178           mySqDist.Append(aVal(anIdx));
 179           myismin.Append(Standard_True);
 180           mypoint.Append(ThePOnC(aParam(anIdx), TheCurveTool::Value(aCurve, aParam(anIdx))));
 181         }
 182         if ((aVal(anIdx) >= aVal(anIdx + 1) &&
 183              aVal(anIdx) >= aVal(anIdx - 1)) ||
 184             (aVal(anIdx) <= aVal(anIdx + 1) &&
 185              aVal(anIdx) <= aVal(anIdx - 1)) )
 186         {
 187           myintuinf = aParam(anIdx - 1);
 188           myintusup = aParam(anIdx + 1);
 189
 190           IntervalPerform(P);
 191         }
 192       }
 193
 194       // Solve on first and last interval.
 195       if (mydist1 > Precision::SquareConfusion())
 196       {
 197         ThePoint aP1, aP2;
 198         TheVector aV1, aV2;
 199         TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower()),     aP1, aV1);
 200         TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower() + 1), aP2, aV2);
 201         TheVector aBase1(P, aP1), aBase2(P, aP2);
 202         Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
 203         Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2
 204
 205         // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
 206         // Necessary condition - when point lies on curve.
 207         if(aVal1 * aVal2 <= 0.0 ||
 208            aBase1.Dot(aBase2) <= 0.0)
 209         {
 210           myintuinf = aParam(aVal.Lower());
 211           myintusup = aParam(aVal.Lower() + 1);
 212           IntervalPerform(P);
 213         }
 214       }
 215
 216       if (mydist2 > Precision::SquareConfusion())
 217       {
 218         ThePoint aP1, aP2;
 219         TheVector aV1, aV2;
 220         TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper() - 1), aP1, aV1);
 221         TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper()),     aP2, aV2);
 222         TheVector aBase1(P, aP1), aBase2(P, aP2);
 223         Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
 224         Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2
 225
 226         // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
 227         // Necessary condition - when point lies on curve.
 228         if(aVal1 * aVal2 <= 0.0 ||
 229            aBase1.Dot(aBase2) <= 0.0)
 230         {
 231           myintuinf = aParam(aVal.Upper() - 1);
 232           myintusup = aParam(aVal.Upper());
 233           IntervalPerform(P);
 234         }
 235       }
 236
 237       mydone = Standard_True;
 238       break;
 239     }
 240   case GeomAbs_OtherCurve:
 241     {
 242       Standard_Boolean IntExtIsDone = Standard_False;
 243       Standard_Boolean IntIsNotValid;
 244       n = TheCurveTool::NbIntervals(aCurve, GeomAbs_C2);
 245       TColStd_Array1OfReal theInter(1, n+1);
 246       Standard_Boolean isPeriodic = TheCurveTool::IsPeriodic(aCurve);
 247       TheCurveTool::Intervals(aCurve, theInter, GeomAbs_C2);
 248       mysample = Max(mysample/n, 17);
 249       TheVector V1;
 250       ThePoint PP;
 251       Standard_Real s1 = 0.0 ;
 252       Standard_Real s2 = 0.0;
 253       myExtPC.Initialize(aCurve);
 254       for (i = 1; i <= n; i++)
 255       {
 256         myintuinf = theInter(i);
 257         myintusup = theInter(i+1);
 258
 259         Standard_Real anInfToCheck = myintuinf;
 260         Standard_Real aSupToCheck = myintusup;
 261
 262         if (isPeriodic) {
 263           Standard_Real aPeriod = TheCurveTool::Period(aCurve);
 264           anInfToCheck = ElCLib::InPeriod(myintuinf, myuinf, myuinf+aPeriod);
 265           aSupToCheck = myintusup+(anInfToCheck-myintuinf);
 266         }
 267         IntIsNotValid = (myuinf > aSupToCheck) || (myusup < anInfToCheck);
 268
 269         if(IntIsNotValid) continue;
 270
 271         if (myuinf >= anInfToCheck) anInfToCheck = myuinf;
 272         if (myusup <= aSupToCheck) aSupToCheck = myusup;
 273         if((aSupToCheck - anInfToCheck) <= mytolu) continue;
 274
 275         if (i != 1)
 276           {
 277           TheCurveTool::D1(aCurve, myintuinf, PP, V1);
 278           s1 = (TheVector(P, PP))*V1;
 279           if (s1*s2 < 0.0) {
 280             mySqDist.Append(PP.SquareDistance(P));
 281             myismin.Append((s1 < 0.0));
 282             mypoint.Append(ThePOnC(myintuinf, PP));
 283           }
 284         }
 285         if (i != n) {
 286           TheCurveTool::D1(aCurve, myintusup, PP, V1);
 287           s2 = (TheVector(P, PP))*V1;
 288         }
 289
 290         IntervalPerform(P);
 291         IntExtIsDone = IntExtIsDone || mydone;
 292       }
 293
 294       mydone = IntExtIsDone;
 295       break;
 296     }
 297   }
 298
 299   // Postprocessing.
 300   if (type == GeomAbs_BSplineCurve ||
 301       type == GeomAbs_OtherCurve)
 302   {
 303     // Additional checking if the point is on the first or last point of the curve
 304     // and does not added yet.
 305     if (mydist1 < Precision::SquareConfusion() ||
 306         mydist2 < Precision::SquareConfusion())
 307     {
 308       Standard_Boolean isFirstAdded = Standard_False;
 309       Standard_Boolean isLastAdded  = Standard_False;
 310       Standard_Integer aNbPoints = mypoint.Length();
 311       for (i = 1; i <= aNbPoints; i++)
 312       {
 313         U = mypoint.Value(i).Parameter();
 314         if (Abs(U - myuinf) < mytolu)
 315           isFirstAdded = Standard_True;
 316         else if (Abs(myusup - U) < mytolu)
 317           isLastAdded = Standard_True;
 318       }
 319       if (!isFirstAdded && mydist1 < Precision::SquareConfusion())
 320       {
 321         mySqDist.Prepend(mydist1);
 322         myismin.Prepend(Standard_True);
 323         mypoint.Prepend(ThePOnC(myuinf, Pf));
 324       }
 325       if (!isLastAdded && mydist2 < Precision::SquareConfusion())
 326       {
 327         mySqDist.Append(mydist2);
 328         myismin.Append(Standard_True);
 329         mypoint.Append(ThePOnC(myusup, Pl));
 330       }
 331       mydone = Standard_True;
 332     }
 333   }
 334   else
 335   {
 336     // In analytical case
 337     mydone = myExtPElC.IsDone();
 338     if (mydone)
 339     {
 340       NbExt = myExtPElC.NbExt();
 341       for (i = 1; i <= NbExt; i++)
 342       {
 343         // Verification de la validite des parametres:
 344         ThePOnC PC = myExtPElC.Point(i);
 345         U = PC.Parameter();
 346         if (TheCurveTool::IsPeriodic(aCurve))
 347         {
 348           U = ElCLib::InPeriod(U, myuinf, myuinf+TheCurveTool::Period(aCurve));
 349         }
 350         if ((U >= myuinf-mytolu) && (U <= myusup+mytolu))
 351         {
 352           PC.SetValues(U, myExtPElC.Point(i).Value());
 353           mySqDist.Append(myExtPElC.SquareDistance(i));
 354           myismin.Append(myExtPElC.IsMin(i));
 355           mypoint.Append(PC);
 356         }
 357       }
 358     }
 359   }
 360 }
 361
 362
 363 //=======================================================================
 364 //function : Initialize
 365 //purpose  :
 366 //=======================================================================
 367
 368 void Extrema_GExtPC::Initialize(const TheCurve&     C,
 369                                 const Standard_Real Uinf,
 370                                 const Standard_Real Usup,
 371                                 const Standard_Real TolF)
 372 {
 373   myC = (Standard_Address)&C;
 374   myintuinf = myuinf = Uinf;
 375   myintusup = myusup = Usup;
 376   mytolf = TolF;
 377   mytolu = TheCurveTool::Resolution(*((TheCurve*)myC), Precision::Confusion());
 378   type = TheCurveTool::GetType(C);
 379   mydone = Standard_False;
 380   mydist1 = RealLast();
 381   mydist2 = RealLast();
 382   mysample = 17;
 383 }
 384
 385
 386 //=======================================================================
 387 //function : IntervalPerform
 388 //purpose  :
 389 //=======================================================================
 390
 391 void Extrema_GExtPC::IntervalPerform(const ThePoint& P)
 392 {
 393   Standard_Integer i;
 394   Standard_Real U;
 395   myExtPC.Initialize(mysample, myintuinf, myintusup, mytolu, mytolf);
 396   myExtPC.Perform(P);
 397   mydone = myExtPC.IsDone();
 398   if (mydone)
 399   {
 400     Standard_Integer NbExt = myExtPC.NbExt();
 401     for (i = 1; i <= NbExt; i++)
 402     {
 403       // Verification de la validite des parametres pour le cas trimme:
 404       ThePOnC PC = myExtPC.Point(i);
 405       U = PC.Parameter();
 406       if (TheCurveTool::IsPeriodic(*((TheCurve*)myC)))
 407       {
 408         U = ElCLib::InPeriod(U, myuinf, myuinf+TheCurveTool::Period(*((TheCurve*)myC)));
 409       }
 410       if ((U >= myuinf - mytolu) && (U <= myusup + mytolu))
 411       {
 412         PC.SetValues(U, PC.Value());
 413         mySqDist.Append(myExtPC.SquareDistance(i));
 414         myismin.Append(myExtPC.IsMin(i));
 415         mypoint.Append(PC);
 416       }
 417     }
 418   }
 419 }
 420
 421
 422
 423
 424 //=======================================================================
 425 //function : Extrema_GExtPC
 426 //purpose  :
 427 //=======================================================================
 428
 429 Extrema_GExtPC::Extrema_GExtPC()
 430 {
 431   myC = 0;
 432   mydone = Standard_False;
 433   mydist1 = RealLast();
 434   mydist2 = RealLast();
 435   mytolu = 0.0;
 436   mytolf = 0.0;
 437   mysample = 17;
 438   myintuinf = myintusup = myuinf = myusup = Precision::Infinite();
 439   type = GeomAbs_OtherCurve;
 440 }
 441
 442 //=======================================================================
 443 //function : Extrema_GExtPC
 444 //purpose  :
 445 //=======================================================================
 446
 447 Extrema_GExtPC::Extrema_GExtPC(const ThePoint&           P,
 448                                const TheCurve&           C,
 449                                const Standard_Real       Uinf,
 450                                const Standard_Real       Usup,
 451                                const Standard_Real       TolF)
 452 {
 453   Initialize(C, Uinf, Usup, TolF);
 454   Perform(P);
 455 }
 456
 457 //=======================================================================
 458 //function : Extrema_GExtPC
 459 //purpose  :
 460 //=======================================================================
 461
 462 Extrema_GExtPC::Extrema_GExtPC(const ThePoint&     P,
 463                                const TheCurve&     C,
 464                                const Standard_Real TolF)
 465 {
 466   Initialize(C, TheCurveTool::FirstParameter(C),
 467              TheCurveTool::LastParameter(C), TolF);
 468   Perform(P);
 469 }
 470
 471
 472 //=======================================================================
 473 //function : IsDone
 474 //purpose  :
 475 //=======================================================================
 476
 477 Standard_Boolean Extrema_GExtPC::IsDone() const
 478 {
 479   return mydone;
 480 }
 481
 482
 483 //=======================================================================
 484 //function : Value
 485 //purpose  :
 486 //=======================================================================
 487
 488 Standard_Real Extrema_GExtPC::SquareDistance(const Standard_Integer N) const
 489 {
 490   if(!mydone) StdFail_NotDone::Raise();
 491   if ((N < 1) || (N > mySqDist.Length())) Standard_OutOfRange::Raise();
 492   return mySqDist.Value(N);
 493 }
 494
 495
 496 //=======================================================================
 497 //function : NbExt
 498 //purpose  :
 499 //=======================================================================
 500
 501 Standard_Integer Extrema_GExtPC::NbExt() const
 502 {
 503   if(!mydone) StdFail_NotDone::Raise();
 504   return mySqDist.Length();
 505 }
 506
 507
 508 //=======================================================================
 509 //function : IsMin
 510 //purpose  :
 511 //=======================================================================
 512
 513 Standard_Boolean Extrema_GExtPC::IsMin(const Standard_Integer N) const
 514 {
 515   if(!mydone) StdFail_NotDone::Raise();
 516   if ((N < 1) || (N > mySqDist.Length())) Standard_OutOfRange::Raise();
 517   return myismin.Value(N);
 518 }
 519
 520
 521
 522 //=======================================================================
 523 //function : Point
 524 //purpose  :
 525 //=======================================================================
 526
 527 const ThePOnC & Extrema_GExtPC::Point(const Standard_Integer N) const
 528 {
 529   if(!mydone) StdFail_NotDone::Raise();
 530   if ((N < 1) || (N > mySqDist.Length())) Standard_OutOfRange::Raise();
 531   return mypoint.Value(N);
 532 }
 533
 534
 535 //=======================================================================
 536 //function : TrimmedDistances
 537 //purpose  :
 538 //=======================================================================
 539
 540 void Extrema_GExtPC::TrimmedSquareDistances(Standard_Real& dist1,
 541                                       Standard_Real& dist2,
 542                                       ThePoint& P1,
 543                                       ThePoint& P2) const
 544 {
 545   dist1 = mydist1;
 546   dist2 = mydist2;
 547   P1 = Pf;
 548   P2 = Pl;
 549 }