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   Standard_Integer NbPoints = (Standard_Integer)(aDiff / Du);
 223   NbPoints = Max(NbPoints, minNbPnts - 1);
 224   Du       = aDiff / NbPoints;
 225
 226   gp_Pnt P;
 227   Standard_Real U = firstu;
 228   for (Standard_Integer i = 1; i <= NbPoints; i++)
 229   {
 230     D0 (C, U, P);
 231     parameters.Append (U);
 232     points    .Append (P);
 233     U += Du;
 234   }
 235   D0 (C, lastu, P);
 236   parameters.Append (lastu);
 237   points    .Append (P);
 238 }
 239
 240
 241 //=======================================================================
 242 //function : PerformCurve
 243 //purpose  : On respecte ll'angle et la fleche, on peut imposer un nombre
 244 //           minimum de points sur un element lineaire
 245 //=======================================================================
 246 void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
 247
 248 {
 249   Standard_Integer i, j;
 250   gp_XYZ V1, V2;
 251   gp_Pnt MiddlePoint, CurrentPoint, LastPoint;
 252   Standard_Real Du, Dusave, MiddleU, L1, L2;
 253
 254   Standard_Real U1       = firstu;
 255   Standard_Real LTol     = Precision::Confusion ();  //protection longueur nulle
 256   Standard_Real ATol     = Precision::Angular ();    //protection angle nul
 257
 258   D0 (C, lastu, LastPoint);
 259
 260   //Initialization du calcul
 261
 262   Standard_Boolean NotDone = Standard_True;
 263   Dusave = (lastu - firstu)*Us3;
 264   Du     = Dusave;
 265   EvaluateDu (C, U1, CurrentPoint, Du, NotDone);
 266   parameters.Append (U1);
 267   points    .Append (CurrentPoint);
 268
 269   // Used to detect "isLine" current bspline and in Du computation in general handling.
 270   Standard_Integer NbInterv = C.NbIntervals(GeomAbs_CN);
 271   TColStd_Array1OfReal Intervs(1, NbInterv+1);
 272   C.Intervals(Intervs, GeomAbs_CN);
 273
 274   if (NotDone) {
 275     //C'est soit une droite, soit une singularite :
 276     V1 = (LastPoint.XYZ() - CurrentPoint.XYZ());
 277     L1 = V1.Modulus ();
 278     if (L1 > LTol)
 279     {
 280       //Si c'est une droite on verifie en calculant minNbPoints :
 281       Standard_Boolean IsLine   = Standard_True;
 282       Standard_Integer NbPoints = (minNbPnts > 3) ? minNbPnts : 3;
 283       ////
 284       Standard_Real param = 0.;
 285       for (i = 1; i <= NbInterv && IsLine; ++i)
 286       {
 287         // Avoid usage intervals out of [firstu, lastu].
 288         if ((Intervs(i+1) < firstu) ||
 289             (Intervs(i)   > lastu))
 290         {
 291           continue;
 292         }
 293         // Fix border points in applicable intervals, to avoid be out of target interval.
 294         if ((Intervs(i)   < firstu) &&
 295             (Intervs(i+1) > firstu))
 296         {
 297           Intervs(i) = firstu;
 298         }
 299         if ((Intervs(i)   < lastu) &&
 300             (Intervs(i+1) > lastu))
 301         {
 302           Intervs(i + 1) = lastu;
 303         }
 304
 305         Standard_Real delta = (Intervs(i+1) - Intervs(i))/NbPoints;
 306         for (j = 1; j <= NbPoints && IsLine; ++j)
 307         {
 308           param = Intervs(i) + j*delta;
 309           D0 (C, param, MiddlePoint);
 310           V2 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
 311           L2 = V2.Modulus ();
 312           if (L2 > LTol)
 313           {
 314             const Standard_Real aAngle = V2.CrossMagnitude(V1)/(L1*L2);
 315             IsLine = (aAngle < ATol);
 316           }
 317         }
 318       }
 319
 320       if (IsLine)
 321       {
 322         parameters.Clear();
 323         points    .Clear();
 324
 325         PerformLinear(C);
 326         return;
 327       }
 328       else
 329       {
 330         //c'etait une singularite on continue :
 331         //Du = Dusave;
 332         EvaluateDu (C, param, MiddlePoint, Du, NotDone);
 333       }
 334     }
 335     else
 336     {
 337       Du = (lastu-firstu)/2.1;
 338       MiddleU = firstu + Du;
 339       D0 (C, MiddleU, MiddlePoint);
 340       V1 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
 341       L1 = V1.Modulus ();
 342       if (L1 < LTol)
 343       {
 344         // L1 < LTol C'est une courbe de longueur nulle, calcul termine :
 345         // on renvoi un segment de 2 points   (protection)
 346         parameters.Append (lastu);
 347         points    .Append (LastPoint);
 348         return;
 349       }
 350     }
 351   }
 352
 353   if (Du > Dusave) Du = Dusave;
 354   else             Dusave = Du;
 355
 356   if (Du < uTol) {
 357     Du = lastu - firstu;
 358     if (Du < uTol) {
 359       parameters.Append (lastu);
 360       points    .Append (LastPoint);
 361       return;
 362     }
 363   }
 364
 365   //Traitement normal pour une courbe
 366   Standard_Boolean MorePoints = Standard_True;
 367   Standard_Real U2            = firstu;
 368   Standard_Real AngleMax      = angularDeflection * 0.5;  //car on prend le point milieu
 369   Standard_Integer aIdx[2] = {Intervs.Lower(), Intervs.Lower()}; // Indexes of intervals of U1 and U2, used to handle non-uniform case.
 370   Standard_Boolean isNeedToCheck = Standard_False;
 371   gp_Pnt aPrevPoint = points.Last();
 372
 373   while (MorePoints) {
 374     aIdx[0] = getIntervalIdx(U1, Intervs, aIdx[0]);
 375     U2 += Du;
 376
 377     if (U2 >= lastu) {                       //Bout de courbe
 378       U2 = lastu;
 379       CurrentPoint = LastPoint;
 380       Du = U2-U1;
 381       Dusave = Du;
 382     }
 383     else D0 (C, U2, CurrentPoint);           //Point suivant
 384
 385     Standard_Real Coef = 0.0, ACoef = 0., FCoef = 0.;
 386     Standard_Boolean Correction, TooLarge, TooSmall;
 387     TooLarge   = Standard_False;
 388     Correction = Standard_True;
 389     TooSmall = Standard_False;
 390
 391     while (Correction) {                     //Ajustement Du
 392       if (isNeedToCheck)
 393       {
 394         aIdx[1] = getIntervalIdx(U2, Intervs, aIdx[0]);
 395         if (aIdx[1] > aIdx[0]) // Jump to another polynom.
 396         {
 397           if (Du > (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3) // Set Du to the smallest value and check deflection on it.
 398           {
 399             Du = (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3;
 400             U2 = U1 + Du;
 401             if (U2 > lastu)
 402               U2 = lastu;
 403             D0 (C, U2, CurrentPoint);
 404           }
 405         }
 406       }
 407       MiddleU = (U1+U2)*0.5;                 //Verif / au point milieu
 408       D0 (C, MiddleU, MiddlePoint);
 409
 410       V1 = (CurrentPoint.XYZ() - aPrevPoint.XYZ()); //Critere de fleche
 411       V2 = (MiddlePoint.XYZ()  - aPrevPoint.XYZ());
 412       L1 = V1.Modulus ();
 413
 414       FCoef = (L1 > myMinLen) ?
 415         V1.CrossMagnitude(V2)/(L1*curvatureDeflection) : 0.0;
 416
 417       V1 = (CurrentPoint.XYZ() - MiddlePoint.XYZ()); //Critere d'angle
 418       L1 = V1.Modulus ();
 419       L2 = V2.Modulus ();
 420       if (L1 > myMinLen && L2 > myMinLen)
 421       {
 422         Standard_Real angg = V1.CrossMagnitude(V2) / (L1 * L2);
 423         ACoef = angg / AngleMax;
 424       }
 425       else
 426         ACoef = 0.0;
 427
 428       //On retient le plus penalisant
 429       Coef = Max(ACoef, FCoef);
 430
 431       if (isNeedToCheck && Coef < 0.55)
 432       {
 433         isNeedToCheck = Standard_False;
 434         Du = Dusave;
 435         U2 = U1 + Du;
 436         if (U2 > lastu)
 437           U2 = lastu;
 438         D0 (C, U2, CurrentPoint);
 439         continue;
 440       }
 441
 442       if (Coef <= 1.0) {
 443         if (Abs (lastu-U2) < uTol) {
 444           parameters.Append (lastu);
 445           points    .Append (LastPoint);
 446           MorePoints = Standard_False;
 447           Correction = Standard_False;
 448         }
 449         else {
 450           if (Coef >= 0.55 || TooLarge) {
 451             parameters.Append (U2);
 452             points    .Append (CurrentPoint);
 453             aPrevPoint = CurrentPoint;
 454             Correction = Standard_False;
 455             isNeedToCheck = Standard_True;
 456           }
 457           else if (TooSmall) {
 458             Correction = Standard_False;
 459             aPrevPoint = CurrentPoint;
 460           }
 461           else {
 462             TooSmall = Standard_True;
 463             //Standard_Real UUU2 = U2;
 464             Du += Min((U2-U1)*(1.-Coef), Du*Us3);
 465
 466             U2 = U1 + Du;
 467             if (U2 > lastu)
 468               U2 = lastu;
 469             D0 (C, U2, CurrentPoint);
 470           }
 471         }
 472       }
 473       else {
 474
 475         if (Coef >= 1.5) {
 476           U2 = MiddleU;
 477           Du  = U2-U1;
 478           CurrentPoint = MiddlePoint;
 479         }
 480         else {
 481           Du*=0.9;
 482           U2 = U1 + Du;
 483           D0 (C, U2, CurrentPoint);
 484           TooLarge = Standard_True;
 485         }
 486
 487       }
 488     }
 489
 490     Du  = U2-U1;
 491
 492     if (MorePoints) {
 493       if (U1 > firstu) {
 494         if (FCoef > ACoef) {
 495           //La fleche est critere de decoupage
 496           EvaluateDu (C, U2, CurrentPoint, Du, NotDone);
 497           if (NotDone) {
 498             Du += (Du-Dusave)*(Du/Dusave);
 499             if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
 500             if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
 501           }
 502         }
 503         else {
 504           //L'angle est le critere de decoupage
 505           Du += (Du-Dusave)*(Du/Dusave);
 506           if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
 507           if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
 508         }
 509       }
 510
 511       if (Du < uTol) {
 512         Du = lastu - U2;
 513         if (Du < uTol) {
 514           parameters.Append (lastu);
 515           points    .Append (LastPoint);
 516           MorePoints = Standard_False;
 517         }
 518         else if (Du*Us3 > uTol) Du*=Us3;
 519       }
 520       U1 = U2;
 521       Dusave = Du;
 522     }
 523   }
 524   //Recalage avant dernier point :
 525   i = points.Length()-1;
 526 //  Real d = points (i).Distance (points (i+1));
 527 // if (Abs(parameters (i) - parameters (i+1))<= 0.000001 || d < Precision::Confusion()) {
 528 //    cout<<"deux points confondus"<<endl;
 529 //    parameters.Remove (i+1);
 530 //    points.Remove (i+1);
 531 //    i--;
 532 //  }
 533   if (i >= 2) {
 534     MiddleU = parameters (i-1);
 535     MiddleU = (lastu + MiddleU)*0.5;
 536     D0 (C, MiddleU, MiddlePoint);
 537     parameters.SetValue (i, MiddleU);
 538     points    .SetValue (i, MiddlePoint);
 539   }
 540
 541   //-- On rajoute des points aux milieux des segments si le nombre
 542   //-- mini de points n'est pas atteint
 543   //--
 544   Standard_Integer Nbp =  points.Length();
 545   Standard_Integer MinNb= (9*minNbPnts)/10;
 546   //if(MinNb<4) MinNb=4;
 547
 548   //-- if(Nbp <  MinNb) { cout<<"\n*"; } else {  cout<<"\n."; }
 549   while(Nbp < MinNb) {
 550     //-- cout<<" \nGCPnts TangentialDeflection : Ajout de Points ("<<Nbp<<" "<<minNbPnts<<" )"<<endl;
 551     for(i=2; i<=Nbp; i++) {
 552       MiddleU = (parameters.Value(i-1)+parameters.Value(i))*0.5;
 553       D0 (C, MiddleU, MiddlePoint);
 554       parameters.InsertBefore(i,MiddleU);
 555       points.InsertBefore(i,MiddlePoint);
 556       Nbp++;
 557       i++;
 558     }
 559   }
 560 }