src/GCPnts/GCPnts_TangentialDeflection.gxx

   1 // Created on: 1996-11-08
   2 // Created by: Jean Claude VAUTHIER
   3 // Copyright (c) 1996-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 <GCPnts_DeflectionType.hxx>
  18 #include <Standard_ConstructionError.hxx>
  19 #include <Precision.hxx>
  20 #include <gp_XYZ.hxx>
  21 #include <gp_Pnt.hxx>
  22 #include <gp_Vec.hxx>
  23 #include <gp.hxx>
  24
  25 #define Us3 0.3333333333333333333333333333
  26
  27 void GCPnts_TangentialDeflection::EvaluateDu (
  28   const TheCurve& C,
  29   const Standard_Real      U,
  30   gp_Pnt&   P,
  31   Standard_Real&     Du,
  32   Standard_Boolean&  NotDone) const {
  33
  34   gp_Vec T, N;
  35   D2 (C, U, P, T, N);
  36   Standard_Real Lt   = T.Magnitude ();
  37   Standard_Real LTol = Precision::Confusion ();
  38   if (Lt > LTol && N.Magnitude () > LTol) {
  39     Standard_Real Lc = N.CrossMagnitude (T);
  40     Standard_Real Ln = Lc/Lt;
  41     if (Ln > LTol) {
  42       Du = sqrt (8.0 * Max(curvatureDeflection, myMinLen) / Ln);
  43       NotDone = Standard_False;
  44     }
  45   }
  46 }
  47
  48
  49 //=======================================================================
  50 //function : GCPnts_TangentialDeflection
  51 //purpose  :
  52 //=======================================================================
  53
  54 GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
  55  const TheCurve&        C,
  56  const Standard_Real    AngularDeflection,
  57  const Standard_Real    CurvatureDeflection,
  58  const Standard_Integer MinimumOfPoints,
  59  const Standard_Real    UTol,
  60  const Standard_Real    theMinLen)
  61
  62 {
  63   Initialize (C,AngularDeflection,CurvatureDeflection,MinimumOfPoints,UTol,theMinLen);
  64 }
  65
  66
  67 //=======================================================================
  68 //function : GCPnts_TangentialDeflection
  69 //purpose  :
  70 //=======================================================================
  71
  72 GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
  73  const TheCurve&        C,
  74  const Standard_Real    FirstParameter,
  75  const Standard_Real    LastParameter,
  76  const Standard_Real    AngularDeflection,
  77  const Standard_Real    CurvatureDeflection,
  78  const Standard_Integer MinimumOfPoints,
  79  const Standard_Real    UTol,
  80  const Standard_Real    theMinLen)
  81
  82 {
  83   Initialize (C,
  84         FirstParameter,
  85         LastParameter,
  86         AngularDeflection,
  87         CurvatureDeflection,
  88         MinimumOfPoints,
  89         UTol, theMinLen);
  90 }
  91
  92
  93
  94 //=======================================================================
  95 //function : Initialize
  96 //purpose  :
  97 //=======================================================================
  98
  99 void GCPnts_TangentialDeflection::Initialize (
 100  const TheCurve&        C,
 101  const Standard_Real    AngularDeflection,
 102  const Standard_Real    CurvatureDeflection,
 103  const Standard_Integer MinimumOfPoints,
 104  const Standard_Real    UTol,
 105  const Standard_Real    theMinLen)
 106
 107 {
 108   Initialize (C,
 109               C.FirstParameter (),
 110               C.LastParameter (),
 111               AngularDeflection,
 112               CurvatureDeflection,
 113               MinimumOfPoints,
 114               UTol, theMinLen);
 115 }
 116
 117
 118 //=======================================================================
 119 //function : Initialize
 120 //purpose  :
 121 //=======================================================================
 122
 123 void GCPnts_TangentialDeflection::Initialize (
 124  const TheCurve&        C,
 125  const Standard_Real    FirstParameter,
 126  const Standard_Real    LastParameter,
 127  const Standard_Real    AngularDeflection,
 128  const Standard_Real    CurvatureDeflection,
 129  const Standard_Integer MinimumOfPoints,
 130  const Standard_Real    UTol,
 131  const Standard_Real    theMinLen)
 132
 133 {
 134
 135   Standard_ConstructionError_Raise_if (CurvatureDeflection <= Precision::Confusion () || AngularDeflection   <= Precision::Angular (), "GCPnts_TangentialDeflection::Initialize - Zero Deflection")
 136
 137  parameters.Clear ();
 138  points    .Clear ();
 139  if (FirstParameter < LastParameter) {
 140    firstu = FirstParameter;
 141    lastu  = LastParameter;
 142  }
 143  else {
 144    lastu  = FirstParameter;
 145    firstu = LastParameter;
 146  }
 147  uTol                = UTol;
 148  angularDeflection   = AngularDeflection;
 149  curvatureDeflection = CurvatureDeflection;
 150  minNbPnts           = Max (MinimumOfPoints, 2);
 151  myMinLen            = Max(theMinLen, Precision::Confusion());
 152
 153  switch (C.GetType()) {
 154
 155  case GeomAbs_Line:
 156    PerformLinear (C);
 157    break;
 158
 159  case GeomAbs_Circle:
 160    PerformCircular (C);
 161    break;
 162
 163  case GeomAbs_BSplineCurve:
 164    {
 165      Handle_TheBSplineCurve BS = C.BSpline() ;
 166      if (BS->NbPoles() == 2 ) PerformLinear (C);
 167      else                     PerformCurve  (C);
 168      break;
 169    }
 170  case GeomAbs_BezierCurve:
 171    {
 172      Handle_TheBezierCurve  BZ = C.Bezier();
 173      if (BZ->NbPoles() == 2) PerformLinear (C);
 174      else                    PerformCurve  (C);
 175      break;
 176    }
 177  default: PerformCurve (C);
 178
 179  }
 180 }
 181
 182
 183 //=======================================================================
 184 //function : PerformLinear
 185 //purpose  :
 186 //=======================================================================
 187
 188 void GCPnts_TangentialDeflection::PerformLinear (const TheCurve& C) {
 189
 190   gp_Pnt P;
 191   D0 (C, firstu, P);
 192   parameters.Append (firstu);
 193   points    .Append (P);
 194   if (minNbPnts > 2) {
 195     Standard_Real Du = (lastu - firstu) / minNbPnts;
 196     Standard_Real U = firstu + Du;
 197     for (Standard_Integer i = 2; i <= minNbPnts; i++) {
 198       D0 (C, U, P);
 199       parameters.Append (U);
 200       points    .Append (P);
 201       U += Du;
 202     }
 203   }
 204   D0 (C, lastu, P);
 205   parameters.Append (lastu);
 206   points    .Append (P);
 207 }
 208
 209 //=======================================================================
 210 //function : PerformCircular
 211 //purpose  :
 212 //=======================================================================
 213
 214 void GCPnts_TangentialDeflection::PerformCircular (const TheCurve& C)
 215 {
 216   // akm 8/01/02 : check the radius before divide by it
 217   Standard_Real dfR = C.Circle().Radius();
 218   Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep(
 219     dfR, curvatureDeflection, angularDeflection, myMinLen);
 220
 221   const Standard_Real aDiff = lastu - firstu;
 222   // Round up number of points to satisfy curvatureDeflection more precisely
 223   Standard_Integer NbPoints = (Standard_Integer)Ceiling(aDiff / Du);
 224   NbPoints = Max(NbPoints, minNbPnts - 1);
 225   Du       = aDiff / NbPoints;
 226
 227   gp_Pnt P;
 228   Standard_Real U = firstu;
 229   for (Standard_Integer i = 1; i <= NbPoints; i++)
 230   {
 231     D0 (C, U, P);
 232     parameters.Append (U);
 233     points    .Append (P);
 234     U += Du;
 235   }
 236   D0 (C, lastu, P);
 237   parameters.Append (lastu);
 238   points    .Append (P);
 239 }
 240
 241
 242 //=======================================================================
 243 //function : PerformCurve
 244 //purpose  : On respecte ll'angle et la fleche, on peut imposer un nombre
 245 //           minimum de points sur un element lineaire
 246 //=======================================================================
 247 void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
 248
 249 {
 250   Standard_Integer i, j;
 251   gp_XYZ V1, V2;
 252   gp_Pnt MiddlePoint, CurrentPoint, LastPoint;
 253   Standard_Real Du, Dusave, MiddleU, L1, L2;
 254
 255   Standard_Real U1       = firstu;
 256   Standard_Real LTol     = Precision::Confusion ();  //protection longueur nulle
 257   Standard_Real ATol     = Precision::Angular ();    //protection angle nul
 258
 259   D0 (C, lastu, LastPoint);
 260
 261   //Initialization du calcul
 262
 263   Standard_Boolean NotDone = Standard_True;
 264   Dusave = (lastu - firstu)*Us3;
 265   Du     = Dusave;
 266   EvaluateDu (C, U1, CurrentPoint, Du, NotDone);
 267   parameters.Append (U1);
 268   points    .Append (CurrentPoint);
 269
 270   // Used to detect "isLine" current bspline and in Du computation in general handling.
 271   Standard_Integer NbInterv = C.NbIntervals(GeomAbs_CN);
 272   TColStd_Array1OfReal Intervs(1, NbInterv+1);
 273   C.Intervals(Intervs, GeomAbs_CN);
 274
 275   if (NotDone) {
 276     //C'est soit une droite, soit une singularite :
 277     V1 = (LastPoint.XYZ() - CurrentPoint.XYZ());
 278     L1 = V1.Modulus ();
 279     if (L1 > LTol)
 280     {
 281       //Si c'est une droite on verifie en calculant minNbPoints :
 282       Standard_Boolean IsLine   = Standard_True;
 283       Standard_Integer NbPoints = (minNbPnts > 3) ? minNbPnts : 3;
 284       ////
 285       Standard_Real param = 0.;
 286       for (i = 1; i <= NbInterv && IsLine; ++i)
 287       {
 288         // Avoid usage intervals out of [firstu, lastu].
 289         if ((Intervs(i+1) < firstu) ||
 290             (Intervs(i)   > lastu))
 291         {
 292           continue;
 293         }
 294         // Fix border points in applicable intervals, to avoid be out of target interval.
 295         if ((Intervs(i)   < firstu) &&
 296             (Intervs(i+1) > firstu))
 297         {
 298           Intervs(i) = firstu;
 299         }
 300         if ((Intervs(i)   < lastu) &&
 301             (Intervs(i+1) > lastu))
 302         {
 303           Intervs(i + 1) = lastu;
 304         }
 305
 306         Standard_Real delta = (Intervs(i+1) - Intervs(i))/NbPoints;
 307         for (j = 1; j <= NbPoints && IsLine; ++j)
 308         {
 309           param = Intervs(i) + j*delta;
 310           D0 (C, param, MiddlePoint);
 311           V2 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
 312           L2 = V2.Modulus ();
 313           if (L2 > LTol)
 314           {
 315             const Standard_Real aAngle = V2.CrossMagnitude(V1)/(L1*L2);
 316             IsLine = (aAngle < ATol);
 317           }
 318         }
 319       }
 320
 321       if (IsLine)
 322       {
 323         parameters.Clear();
 324         points    .Clear();
 325
 326         PerformLinear(C);
 327         return;
 328       }
 329       else
 330       {
 331         //c'etait une singularite on continue :
 332         //Du = Dusave;
 333         EvaluateDu (C, param, MiddlePoint, Du, NotDone);
 334       }
 335     }
 336     else
 337     {
 338       Du = (lastu-firstu)/2.1;
 339       MiddleU = firstu + Du;
 340       D0 (C, MiddleU, MiddlePoint);
 341       V1 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
 342       L1 = V1.Modulus ();
 343       if (L1 < LTol)
 344       {
 345         // L1 < LTol C'est une courbe de longueur nulle, calcul termine :
 346         // on renvoi un segment de 2 points   (protection)
 347         parameters.Append (lastu);
 348         points    .Append (LastPoint);
 349         return;
 350       }
 351     }
 352   }
 353
 354   if (Du > Dusave) Du = Dusave;
 355   else             Dusave = Du;
 356
 357   if (Du < uTol) {
 358     Du = lastu - firstu;
 359     if (Du < uTol) {
 360       parameters.Append (lastu);
 361       points    .Append (LastPoint);
 362       return;
 363     }
 364   }
 365
 366   //Traitement normal pour une courbe
 367   Standard_Boolean MorePoints = Standard_True;
 368   Standard_Real U2            = firstu;
 369   Standard_Real AngleMax      = angularDeflection * 0.5;  //car on prend le point milieu
 370   Standard_Integer aIdx[2] = {Intervs.Lower(), Intervs.Lower()}; // Indexes of intervals of U1 and U2, used to handle non-uniform case.
 371   Standard_Boolean isNeedToCheck = Standard_False;
 372   gp_Pnt aPrevPoint = points.Last();
 373
 374   while (MorePoints) {
 375     aIdx[0] = getIntervalIdx(U1, Intervs, aIdx[0]);
 376     U2 += Du;
 377
 378     if (U2 >= lastu) {                       //Bout de courbe
 379       U2 = lastu;
 380       CurrentPoint = LastPoint;
 381       Du = U2-U1;
 382       Dusave = Du;
 383     }
 384     else D0 (C, U2, CurrentPoint);           //Point suivant
 385
 386     Standard_Real Coef = 0.0, ACoef = 0., FCoef = 0.;
 387     Standard_Boolean Correction, TooLarge, TooSmall;
 388     TooLarge   = Standard_False;
 389     Correction = Standard_True;
 390     TooSmall = Standard_False;
 391
 392     while (Correction) {                     //Ajustement Du
 393       if (isNeedToCheck)
 394       {
 395         aIdx[1] = getIntervalIdx(U2, Intervs, aIdx[0]);
 396         if (aIdx[1] > aIdx[0]) // Jump to another polynom.
 397         {
 398           if (Du > (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3) // Set Du to the smallest value and check deflection on it.
 399           {
 400             Du = (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3;
 401             U2 = U1 + Du;
 402             if (U2 > lastu)
 403               U2 = lastu;
 404             D0 (C, U2, CurrentPoint);
 405           }
 406         }
 407       }
 408       MiddleU = (U1+U2)*0.5;                 //Verif / au point milieu
 409       D0 (C, MiddleU, MiddlePoint);
 410
 411       V1 = (CurrentPoint.XYZ() - aPrevPoint.XYZ()); //Critere de fleche
 412       V2 = (MiddlePoint.XYZ()  - aPrevPoint.XYZ());
 413       L1 = V1.Modulus ();
 414
 415       FCoef = (L1 > myMinLen) ?
 416         V1.CrossMagnitude(V2)/(L1*curvatureDeflection) : 0.0;
 417
 418       V1 = (CurrentPoint.XYZ() - MiddlePoint.XYZ()); //Critere d'angle
 419       L1 = V1.Modulus ();
 420       L2 = V2.Modulus ();
 421       if (L1 > myMinLen && L2 > myMinLen)
 422       {
 423         Standard_Real angg = V1.CrossMagnitude(V2) / (L1 * L2);
 424         ACoef = angg / AngleMax;
 425       }
 426       else
 427         ACoef = 0.0;
 428
 429       //On retient le plus penalisant
 430       Coef = Max(ACoef, FCoef);
 431
 432       if (isNeedToCheck && Coef < 0.55)
 433       {
 434         isNeedToCheck = Standard_False;
 435         Du = Dusave;
 436         U2 = U1 + Du;
 437         if (U2 > lastu)
 438           U2 = lastu;
 439         D0 (C, U2, CurrentPoint);
 440         continue;
 441       }
 442
 443       if (Coef <= 1.0) {
 444         if (Abs (lastu-U2) < uTol) {
 445           parameters.Append (lastu);
 446           points    .Append (LastPoint);
 447           MorePoints = Standard_False;
 448           Correction = Standard_False;
 449         }
 450         else {
 451           if (Coef >= 0.55 || TooLarge) {
 452             parameters.Append (U2);
 453             points    .Append (CurrentPoint);
 454             aPrevPoint = CurrentPoint;
 455             Correction = Standard_False;
 456             isNeedToCheck = Standard_True;
 457           }
 458           else if (TooSmall) {
 459             Correction = Standard_False;
 460             aPrevPoint = CurrentPoint;
 461           }
 462           else {
 463             TooSmall = Standard_True;
 464             //Standard_Real UUU2 = U2;
 465             Du += Min((U2-U1)*(1.-Coef), Du*Us3);
 466
 467             U2 = U1 + Du;
 468             if (U2 > lastu)
 469               U2 = lastu;
 470             D0 (C, U2, CurrentPoint);
 471           }
 472         }
 473       }
 474       else {
 475
 476         if (Coef >= 1.5) {
 477           if (!aPrevPoint.IsEqual(points.Last(), Precision::Confusion()))
 478           {
 479             parameters.Append (U1);
 480             points    .Append (aPrevPoint);
 481           }
 482           U2 = MiddleU;
 483           Du  = U2-U1;
 484           CurrentPoint = MiddlePoint;
 485         }
 486         else {
 487           Du*=0.9;
 488           U2 = U1 + Du;
 489           D0 (C, U2, CurrentPoint);
 490           TooLarge = Standard_True;
 491         }
 492
 493       }
 494     }
 495
 496     Du  = U2-U1;
 497
 498     if (MorePoints) {
 499       if (U1 > firstu) {
 500         if (FCoef > ACoef) {
 501           //La fleche est critere de decoupage
 502           EvaluateDu (C, U2, CurrentPoint, Du, NotDone);
 503           if (NotDone) {
 504             Du += (Du-Dusave)*(Du/Dusave);
 505             if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
 506             if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
 507           }
 508         }
 509         else {
 510           //L'angle est le critere de decoupage
 511           Du += (Du-Dusave)*(Du/Dusave);
 512           if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
 513           if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
 514         }
 515       }
 516
 517       if (Du < uTol) {
 518         Du = lastu - U2;
 519         if (Du < uTol) {
 520           parameters.Append (lastu);
 521           points    .Append (LastPoint);
 522           MorePoints = Standard_False;
 523         }
 524         else if (Du*Us3 > uTol) Du*=Us3;
 525       }
 526       U1 = U2;
 527       Dusave = Du;
 528     }
 529   }
 530   //Recalage avant dernier point :
 531   i = points.Length()-1;
 532 //  Real d = points (i).Distance (points (i+1));
 533 // if (Abs(parameters (i) - parameters (i+1))<= 0.000001 || d < Precision::Confusion()) {
 534 //    cout<<"deux points confondus"<<endl;
 535 //    parameters.Remove (i+1);
 536 //    points.Remove (i+1);
 537 //    i--;
 538 //  }
 539   if (i >= 2) {
 540     MiddleU = parameters (i-1);
 541     MiddleU = (lastu + MiddleU)*0.5;
 542     D0 (C, MiddleU, MiddlePoint);
 543     parameters.SetValue (i, MiddleU);
 544     points    .SetValue (i, MiddlePoint);
 545   }
 546
 547   //-- On rajoute des points aux milieux des segments si le nombre
 548   //-- mini de points n'est pas atteint
 549   //--
 550   Standard_Integer Nbp =  points.Length();
 551   Standard_Integer MinNb= (9*minNbPnts)/10;
 552   //if(MinNb<4) MinNb=4;
 553
 554   //-- if(Nbp <  MinNb) { cout<<"\n*"; } else {  cout<<"\n."; }
 555   while(Nbp < MinNb) {
 556     //-- cout<<" \nGCPnts TangentialDeflection : Ajout de Points ("<<Nbp<<" "<<minNbPnts<<" )"<<endl;
 557     for(i=2; i<=Nbp; i++) {
 558       MiddleU = (parameters.Value(i-1)+parameters.Value(i))*0.5;
 559       D0 (C, MiddleU, MiddlePoint);
 560       parameters.InsertBefore(i,MiddleU);
 561       points.InsertBefore(i,MiddlePoint);
 562       Nbp++;
 563       i++;
 564     }
 565   }
 566 }