src/GeomEvaluator/GeomEvaluator_OffsetSurface.cxx

   1 // Created on: 2015-09-21
   2 // Copyright (c) 2015 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 #include <GeomEvaluator_OffsetSurface.hxx>
  16
  17 #include <GeomAdaptor_HSurface.hxx>
  18 #include <CSLib.hxx>
  19 #include <CSLib_NormalStatus.hxx>
  20 #include <Geom_BezierSurface.hxx>
  21 #include <Geom_BSplineSurface.hxx>
  22 #include <Geom_UndefinedValue.hxx>
  23 #include <GeomAdaptor_HSurface.hxx>
  24 #include <gp_Vec.hxx>
  25 #include <Standard_RangeError.hxx>
  26 #include <Standard_NumericError.hxx>
  27
  28 IMPLEMENT_STANDARD_RTTIEXT(GeomEvaluator_OffsetSurface,GeomEvaluator_Surface)
  29
  30 namespace {
  31
  32 // tolerance for considering derivative to be null
  33 const Standard_Real the_D1MagTol = 1.e-9;
  34
  35 // If calculation of normal fails, try shifting the point towards the center
  36 // of the parametric space of the surface, in the hope that derivatives
  37 // are better defined there.
  38 //
  39 // This shift is iterative, starting with Precision::PConfusion()
  40 // and increasing by multiple of 2 on each step.
  41 //
  42 // NB: temporarily this is made as static function and not class method,
  43 // hence code duplications
  44 static Standard_Boolean shiftPoint (const Standard_Real theUStart, const Standard_Real theVStart,
  45                                     Standard_Real& theU, Standard_Real& theV,
  46                                     const Handle(Geom_Surface)& theSurf,
  47                                     const Handle(GeomAdaptor_HSurface)& theAdaptor,
  48                                     const gp_Vec& theD1U, const gp_Vec& theD1V)
  49 {
  50   // Get parametric bounds and closure status
  51   Standard_Real aUMin, aUMax, aVMin, aVMax;
  52   Standard_Boolean isUPeriodic, isVPeriodic;
  53   if (! theSurf.IsNull())
  54   {
  55     theSurf->Bounds (aUMin, aUMax, aVMin, aVMax);
  56     isUPeriodic = theSurf->IsUPeriodic();
  57     isVPeriodic = theSurf->IsVPeriodic();
  58   }
  59   else
  60   {
  61     aUMin = theAdaptor->FirstUParameter();
  62     aUMax = theAdaptor->LastUParameter();
  63     aVMin = theAdaptor->FirstVParameter();
  64     aVMax = theAdaptor->LastVParameter();
  65     isUPeriodic = theAdaptor->IsUPeriodic();
  66     isVPeriodic = theAdaptor->IsVPeriodic();
  67   }
  68
  69   // check if either U or V is singular (normally one of them is)
  70   Standard_Boolean isUSingular = (theD1U.SquareMagnitude() < the_D1MagTol * the_D1MagTol);
  71   Standard_Boolean isVSingular = (theD1V.SquareMagnitude() < the_D1MagTol * the_D1MagTol);
  72
  73   // compute vector to shift from start point to center of the surface;
  74   // if surface is periodic or singular in some direction, take shift in that direction zero
  75   Standard_Real aDirU = (isUPeriodic || (isUSingular && ! isVSingular) ? 0. : 0.5 * (aUMin + aUMax) - theUStart);
  76   Standard_Real aDirV = (isVPeriodic || (isVSingular && ! isUSingular) ? 0. : 0.5 * (aVMin + aVMax) - theVStart);
  77   Standard_Real aDist = Sqrt (aDirU * aDirU + aDirV * aDirV);
  78
  79   // shift current point from its current position towards center, by value of twice
  80   // current distance from it to start (but not less than Precision::PConfusion());
  81   // fail if center is overpassed.
  82   Standard_Real aDU = theU - theUStart;
  83   Standard_Real aDV = theV - theVStart;
  84   Standard_Real aStep = Max (2. * Sqrt (aDU * aDU + aDV * aDV), Precision::PConfusion());
  85   if (aStep >= aDist)
  86   {
  87     return Standard_False;
  88   }
  89
  90   aStep /= aDist;
  91   theU += aDirU * aStep;
  92   theV += aDirV * aStep;
  93
  94 //  std::cout << "Normal calculation failed at (" << theUStart << ", " << theVStart << "), shifting to (" << theU << ", " << theV << ")" << std::endl;
  95
  96   return Standard_True;
  97 }
  98
  99 // auxiliary function
 100 template<class SurfOrAdapt>
 101 static void derivatives(Standard_Integer theMaxOrder,
 102                         Standard_Integer theMinOrder,
 103                         const Standard_Real theU,
 104                         const Standard_Real theV,
 105                         const SurfOrAdapt& theBasisSurf,
 106                         const Standard_Integer theNU,
 107                         const Standard_Integer theNV,
 108                         const Standard_Boolean theAlongU,
 109                         const Standard_Boolean theAlongV,
 110                         const Handle(Geom_BSplineSurface)& theL,
 111                         TColgp_Array2OfVec& theDerNUV,
 112                         TColgp_Array2OfVec& theDerSurf)
 113 {
 114   Standard_Integer i, j;
 115   gp_Pnt P;
 116   gp_Vec DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3UUV, DL3UVV, DL3V;
 117
 118   if (theAlongU || theAlongV)
 119   {
 120     theMaxOrder = 0;
 121     TColgp_Array2OfVec DerSurfL(0, theMaxOrder + theNU + 1, 0, theMaxOrder + theNV + 1);
 122     switch (theMinOrder)
 123     {
 124     case 1:
 125       theL->D1(theU, theV, P, DL1U, DL1V);
 126       DerSurfL.SetValue(1, 0, DL1U);
 127       DerSurfL.SetValue(0, 1, DL1V);
 128       break;
 129     case 2:
 130       theL->D2(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV);
 131       DerSurfL.SetValue(1, 0, DL1U);
 132       DerSurfL.SetValue(0, 1, DL1V);
 133       DerSurfL.SetValue(1, 1, DL2UV);
 134       DerSurfL.SetValue(2, 0, DL2U);
 135       DerSurfL.SetValue(0, 2, DL2V);
 136       break;
 137     case 3:
 138       theL->D3(theU, theV, P, DL1U, DL1V, DL2U, DL2V, DL2UV, DL3U, DL3V, DL3UUV, DL3UVV);
 139       DerSurfL.SetValue(1, 0, DL1U);
 140       DerSurfL.SetValue(0, 1, DL1V);
 141       DerSurfL.SetValue(1, 1, DL2UV);
 142       DerSurfL.SetValue(2, 0, DL2U);
 143       DerSurfL.SetValue(0, 2, DL2V);
 144       DerSurfL.SetValue(3, 0, DL3U);
 145       DerSurfL.SetValue(2, 1, DL3UUV);
 146       DerSurfL.SetValue(1, 2, DL3UVV);
 147       DerSurfL.SetValue(0, 3, DL3V);
 148       break;
 149     default:
 150       break;
 151     }
 152
 153     if (theNU <= theNV)
 154     {
 155       for (i = 0; i <= theMaxOrder + 1 + theNU; i++)
 156         for (j = i; j <= theMaxOrder + theNV + 1; j++)
 157           if (i + j > theMinOrder)
 158           {
 159             DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
 160             theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
 161             if (i != j && j <= theNU + 1)
 162             {
 163               theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
 164               DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
 165             }
 166           }
 167     }
 168     else
 169     {
 170       for (j = 0; j <= theMaxOrder + 1 + theNV; j++)
 171         for (i = j; i <= theMaxOrder + theNU + 1; i++)
 172           if (i + j > theMinOrder)
 173           {
 174             DerSurfL.SetValue(i, j, theL->DN(theU, theV, i, j));
 175             theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
 176             if (i != j && i <= theNV + 1)
 177             {
 178               theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
 179               DerSurfL.SetValue(j, i, theL->DN(theU, theV, j, i));
 180             }
 181           }
 182     }
 183     for (i = 0; i <= theMaxOrder + theNU; i++)
 184       for (j = 0; j <= theMaxOrder + theNV; j++)
 185       {
 186         if (theAlongU)
 187           theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, DerSurfL, theDerSurf));
 188         if (theAlongV)
 189           theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf, DerSurfL));
 190       }
 191   }
 192   else
 193   {
 194     for (i = 0; i <= theMaxOrder + theNU+ 1; i++)
 195     {
 196       for (j = i; j <= theMaxOrder + theNV + 1; j++)
 197       {
 198         if (i + j > theMinOrder)
 199         {
 200           theDerSurf.SetValue(i, j, theBasisSurf->DN(theU, theV, i, j));
 201           if (i != j
 202            && j <= theDerSurf.UpperRow()
 203            && i <= theDerSurf.UpperCol())
 204           {
 205             theDerSurf.SetValue(j, i, theBasisSurf->DN(theU, theV, j, i));
 206           }
 207         }
 208       }
 209     }
 210     for (i = 0; i <= theMaxOrder + theNU; i++)
 211       for (j = 0; j <= theMaxOrder + theNV; j++)
 212         theDerNUV.SetValue(i, j, CSLib::DNNUV(i, j, theDerSurf));
 213   }
 214 }
 215
 216 inline Standard_Boolean IsInfiniteCoord (const gp_Vec& theVec)
 217 {
 218   return Precision::IsInfinite (theVec.X()) ||
 219          Precision::IsInfinite (theVec.Y()) ||
 220          Precision::IsInfinite (theVec.Z());
 221 }
 222
 223 inline void CheckInfinite (const gp_Vec& theVecU, const gp_Vec& theVecV)
 224 {
 225   if (IsInfiniteCoord (theVecU) || IsInfiniteCoord (theVecV))
 226   {
 227     throw Standard_NumericError ("GeomEvaluator_OffsetSurface: Evaluation of infinite parameters");
 228   }
 229 }
 230
 231 } // end of anonymous namespace
 232
 233 GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface(
 234         const Handle(Geom_Surface)& theBase,
 235         const Standard_Real theOffset,
 236         const Handle(Geom_OsculatingSurface)& theOscSurf)
 237   : GeomEvaluator_Surface(),
 238     myBaseSurf(theBase),
 239     myOffset(theOffset),
 240     myOscSurf(theOscSurf)
 241 {
 242   if (!myOscSurf.IsNull())
 243     return; // osculating surface already exists
 244
 245   // Create osculating surface for B-spline and Besier surfaces only
 246   if (myBaseSurf->IsKind(STANDARD_TYPE(Geom_BSplineSurface)) ||
 247       myBaseSurf->IsKind(STANDARD_TYPE(Geom_BezierSurface)))
 248     myOscSurf = new Geom_OsculatingSurface(myBaseSurf, Precision::Confusion());
 249 }
 250
 251 GeomEvaluator_OffsetSurface::GeomEvaluator_OffsetSurface(
 252         const Handle(GeomAdaptor_HSurface)& theBase,
 253         const Standard_Real theOffset,
 254         const Handle(Geom_OsculatingSurface)& theOscSurf)
 255   : GeomEvaluator_Surface(),
 256     myBaseAdaptor(theBase),
 257     myOffset(theOffset),
 258     myOscSurf(theOscSurf)
 259 {
 260 }
 261
 262 void GeomEvaluator_OffsetSurface::D0(
 263     const Standard_Real theU, const Standard_Real theV, gp_Pnt& theValue) const
 264 {
 265   Standard_Real aU = theU, aV = theV;
 266   for (;;)
 267   {
 268     gp_Vec aD1U, aD1V;
 269     BaseD1 (aU, aV, theValue, aD1U, aD1V);
 270
 271     CheckInfinite (aD1U, aD1V);
 272
 273     try
 274     {
 275       CalculateD0(aU, aV, theValue, aD1U, aD1V);
 276       break;
 277     }
 278     catch (Geom_UndefinedValue&)
 279     {
 280       // if failed at parametric boundary, try taking derivative at shifted point
 281       if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, aD1U, aD1V))
 282       {
 283         throw;
 284       }
 285     }
 286   }
 287 }
 288
 289 void GeomEvaluator_OffsetSurface::D1(
 290     const Standard_Real theU, const Standard_Real theV,
 291     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V) const
 292 {
 293   Standard_Real aU = theU, aV = theV;
 294   for (;;)
 295   {
 296     gp_Vec aD2U, aD2V, aD2UV;
 297     BaseD2 (aU, aV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
 298
 299     CheckInfinite (theD1U, theD1V);
 300
 301     try
 302     {
 303       CalculateD1(aU, aV, theValue, theD1U, theD1V, aD2U, aD2V, aD2UV);
 304       break;
 305     }
 306     catch (Geom_UndefinedValue&)
 307     {
 308       // if failed at parametric boundary, try taking derivative at shifted point
 309       if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V))
 310       {
 311         throw;
 312       }
 313     }
 314   }
 315 }
 316
 317 void GeomEvaluator_OffsetSurface::D2(
 318     const Standard_Real theU, const Standard_Real theV,
 319     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 320     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV) const
 321 {
 322   Standard_Real aU = theU, aV = theV;
 323   for (;;)
 324   {
 325     gp_Vec aD3U, aD3V, aD3UUV, aD3UVV;
 326     BaseD3 (aU, aV, theValue, theD1U, theD1V,
 327             theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
 328
 329     CheckInfinite (theD1U, theD1V);
 330
 331     try
 332     {
 333       CalculateD2 (aU, aV, theValue, theD1U, theD1V,
 334                    theD2U, theD2V, theD2UV, aD3U, aD3V, aD3UUV, aD3UVV);
 335       break;
 336     }
 337     catch (Geom_UndefinedValue&)
 338     {
 339       // if failed at parametric boundary, try taking derivative at shifted point
 340       if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V))
 341       {
 342         throw;
 343       }
 344     }
 345   }
 346 }
 347
 348 void GeomEvaluator_OffsetSurface::D3(
 349     const Standard_Real theU, const Standard_Real theV,
 350     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 351     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
 352     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
 353 {
 354   Standard_Real aU = theU, aV = theV;
 355   for (;;)
 356   {
 357     BaseD3 (aU, aV, theValue, theD1U, theD1V,
 358             theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
 359
 360     CheckInfinite (theD1U, theD1V);
 361
 362     try
 363     {
 364       CalculateD3 (aU, aV, theValue, theD1U, theD1V,
 365                    theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
 366       break;
 367     }
 368     catch (Geom_UndefinedValue&)
 369     {
 370       // if failed at parametric boundary, try taking derivative at shifted point
 371       if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, theD1U, theD1V))
 372       {
 373         throw;
 374       }
 375     }
 376   }
 377 }
 378
 379 gp_Vec GeomEvaluator_OffsetSurface::DN(
 380     const Standard_Real theU, const Standard_Real theV,
 381     const Standard_Integer theDerU, const Standard_Integer theDerV) const
 382 {
 383   Standard_RangeError_Raise_if(theDerU < 0, "GeomEvaluator_OffsetSurface::DN(): theDerU < 0");
 384   Standard_RangeError_Raise_if(theDerV < 0, "GeomEvaluator_OffsetSurface::DN(): theDerV < 0");
 385   Standard_RangeError_Raise_if(theDerU + theDerV < 1,
 386       "GeomEvaluator_OffsetSurface::DN(): theDerU + theDerV < 1");
 387
 388   Standard_Real aU = theU, aV = theV;
 389   for (;;)
 390   {
 391     gp_Pnt aP;
 392     gp_Vec aD1U, aD1V;
 393     BaseD1 (aU, aV, aP, aD1U, aD1V);
 394
 395     CheckInfinite (aD1U, aD1V);
 396
 397     try
 398     {
 399       return CalculateDN (aU, aV, theDerU, theDerV, aD1U, aD1V);
 400     }
 401     catch (Geom_UndefinedValue&)
 402     {
 403       // if failed at parametric boundary, try taking derivative at shifted point
 404       if (! shiftPoint (theU, theV, aU, aV, myBaseSurf, myBaseAdaptor, aD1U, aD1V))
 405       {
 406         throw;
 407       }
 408     }
 409   }
 410 }
 411
 412
 413 void GeomEvaluator_OffsetSurface::BaseD0(const Standard_Real theU, const Standard_Real theV,
 414                                          gp_Pnt& theValue) const
 415 {
 416   if (!myBaseAdaptor.IsNull())
 417     myBaseAdaptor->D0(theU, theV, theValue);
 418   else
 419     myBaseSurf->D0(theU, theV, theValue);
 420 }
 421
 422 void GeomEvaluator_OffsetSurface::BaseD1(const Standard_Real theU, const Standard_Real theV,
 423                                          gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V) const
 424 {
 425   if (!myBaseAdaptor.IsNull())
 426     myBaseAdaptor->D1(theU, theV, theValue, theD1U, theD1V);
 427   else
 428     myBaseSurf->D1(theU, theV, theValue, theD1U, theD1V);
 429 }
 430
 431 void GeomEvaluator_OffsetSurface::BaseD2(const Standard_Real theU, const Standard_Real theV,
 432                                          gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 433                                          gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV) const
 434 {
 435   if (!myBaseAdaptor.IsNull())
 436     myBaseAdaptor->D2(theU, theV, theValue, theD1U, theD1V, theD2U, theD2V, theD2UV);
 437   else
 438     myBaseSurf->D2(theU, theV, theValue, theD1U, theD1V, theD2U, theD2V, theD2UV);
 439 }
 440
 441 void GeomEvaluator_OffsetSurface::BaseD3(
 442     const Standard_Real theU, const Standard_Real theV,
 443     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 444     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
 445     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
 446 {
 447   if (!myBaseAdaptor.IsNull())
 448     myBaseAdaptor->D3(theU, theV, theValue, theD1U, theD1V,
 449     theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
 450   else
 451     myBaseSurf->D3(theU, theV, theValue, theD1U, theD1V,
 452     theD2U, theD2V, theD2UV, theD3U, theD3V, theD3UUV, theD3UVV);
 453 }
 454
 455
 456 void GeomEvaluator_OffsetSurface::CalculateD0(
 457     const Standard_Real theU, const Standard_Real theV,
 458     gp_Pnt& theValue,
 459     const gp_Vec& theD1U, const gp_Vec& theD1V) const
 460 {
 461   // Normalize derivatives before normal calculation because it gives more stable result.
 462   // There will be normalized only derivatives greater than 1.0 to avoid differences in last significant digit
 463   gp_Vec aD1U(theD1U);
 464   gp_Vec aD1V(theD1V);
 465   Standard_Real aD1UNorm2 = aD1U.SquareMagnitude();
 466   Standard_Real aD1VNorm2 = aD1V.SquareMagnitude();
 467   if (aD1UNorm2 > 1.0)
 468     aD1U /= Sqrt(aD1UNorm2);
 469   if (aD1VNorm2 > 1.0)
 470     aD1V /= Sqrt(aD1VNorm2);
 471
 472   gp_Vec aNorm = aD1U.Crossed(aD1V);
 473   if (aNorm.SquareMagnitude() > the_D1MagTol * the_D1MagTol)
 474   {
 475     // Non singular case. Simple computations.
 476     aNorm.Normalize();
 477     theValue.SetXYZ(theValue.XYZ() + myOffset * aNorm.XYZ());
 478   }
 479   else
 480   {
 481     const Standard_Integer MaxOrder = 3;
 482
 483     Handle(Geom_BSplineSurface) L;
 484     Standard_Boolean isOpposite = Standard_False;
 485     Standard_Boolean AlongU = Standard_False;
 486     Standard_Boolean AlongV = Standard_False;
 487     if (!myOscSurf.IsNull())
 488     {
 489       AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
 490       AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
 491     }
 492     const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
 493
 494     TColgp_Array2OfVec DerNUV(0, MaxOrder, 0, MaxOrder);
 495     TColgp_Array2OfVec DerSurf(0, MaxOrder + 1, 0, MaxOrder + 1);
 496     Standard_Integer OrderU, OrderV;
 497     Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
 498     Bounds(Umin, Umax, Vmin, Vmax);
 499
 500     DerSurf.SetValue(1, 0, theD1U);
 501     DerSurf.SetValue(0, 1, theD1V);
 502     if (!myBaseSurf.IsNull())
 503       derivatives(MaxOrder, 1, theU, theV, myBaseSurf, 0, 0, AlongU, AlongV, L, DerNUV, DerSurf);
 504     else
 505       derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, 0, 0, AlongU, AlongV, L, DerNUV, DerSurf);
 506
 507     gp_Dir Normal;
 508     CSLib_NormalStatus NStatus = CSLib_Singular;
 509     CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
 510                   NStatus, Normal, OrderU, OrderV);
 511     if (NStatus == CSLib_InfinityOfSolutions)
 512     {
 513       // Replace zero derivative and try to calculate normal
 514       gp_Vec aNewDU = theD1U;
 515       gp_Vec aNewDV = theD1V;
 516       if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, the_D1MagTol * the_D1MagTol))
 517         CSLib::Normal(aNewDU, aNewDV, the_D1MagTol, NStatus, Normal);
 518     }
 519
 520     if (NStatus != CSLib_Defined)
 521       throw Geom_UndefinedValue(
 522           "GeomEvaluator_OffsetSurface::CalculateD0(): Unable to calculate normal");
 523
 524     theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
 525   }
 526 }
 527
 528 void GeomEvaluator_OffsetSurface::CalculateD1(
 529     const Standard_Real theU, const Standard_Real theV,
 530     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 531     const gp_Vec& theD2U, const gp_Vec& theD2V, const gp_Vec& theD2UV) const
 532 {
 533   // Check offset side.
 534   Handle(Geom_BSplineSurface) L;
 535   Standard_Boolean isOpposite = Standard_False;
 536   Standard_Boolean AlongU = Standard_False;
 537   Standard_Boolean AlongV = Standard_False;
 538
 539   // Normalize derivatives before normal calculation because it gives more stable result.
 540   // There will be normalized only derivatives greater than 1.0 to avoid differences in last significant digit
 541   gp_Vec aD1U(theD1U);
 542   gp_Vec aD1V(theD1V);
 543   Standard_Real aD1UNorm2 = aD1U.SquareMagnitude();
 544   Standard_Real aD1VNorm2 = aD1V.SquareMagnitude();
 545   if (aD1UNorm2 > 1.0)
 546     aD1U /= Sqrt(aD1UNorm2);
 547   if (aD1VNorm2 > 1.0)
 548     aD1V /= Sqrt(aD1VNorm2);
 549
 550   Standard_Boolean isSingular = Standard_False;
 551   Standard_Integer MaxOrder = 0;
 552   gp_Vec aNorm = aD1U.Crossed(aD1V);
 553   if (aNorm.SquareMagnitude() <= the_D1MagTol * the_D1MagTol)
 554   {
 555     if (!myOscSurf.IsNull())
 556     {
 557       AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
 558       AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
 559     }
 560
 561     MaxOrder = 3;
 562     isSingular = Standard_True;
 563   }
 564
 565   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
 566
 567   if (!isSingular)
 568   {
 569     // AlongU or AlongV leads to more complex D1 computation
 570     // Try to compute D0 and D1 much simpler
 571     aNorm.Normalize();
 572     theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * aNorm.XYZ());
 573
 574     gp_Vec aN0(aNorm.XYZ()), aN1U, aN1V;
 575     Standard_Real aScale = (theD1U^theD1V).Dot(aN0);
 576     aN1U.SetX(theD2U.Y() * theD1V.Z() + theD1U.Y() * theD2UV.Z()
 577       - theD2U.Z() * theD1V.Y() - theD1U.Z() * theD2UV.Y());
 578     aN1U.SetY((theD2U.X() * theD1V.Z() + theD1U.X() * theD2UV.Z()
 579       - theD2U.Z() * theD1V.X() - theD1U.Z() * theD2UV.X()) * -1.0);
 580     aN1U.SetZ(theD2U.X() * theD1V.Y() + theD1U.X() * theD2UV.Y()
 581       - theD2U.Y() * theD1V.X() - theD1U.Y() * theD2UV.X());
 582     Standard_Real aScaleU = aN1U.Dot(aN0);
 583     aN1U.Subtract(aScaleU * aN0);
 584     aN1U /= aScale;
 585
 586     aN1V.SetX(theD2UV.Y() * theD1V.Z() + theD2V.Z() * theD1U.Y()
 587       - theD2UV.Z() * theD1V.Y() - theD2V.Y() * theD1U.Z());
 588     aN1V.SetY((theD2UV.X() * theD1V.Z() + theD2V.Z() * theD1U.X()
 589       - theD2UV.Z() * theD1V.X() - theD2V.X() * theD1U.Z()) * -1.0);
 590     aN1V.SetZ(theD2UV.X() * theD1V.Y() + theD2V.Y() * theD1U.X()
 591       - theD2UV.Y() * theD1V.X() - theD2V.X() * theD1U.Y());
 592     Standard_Real aScaleV = aN1V.Dot(aN0);
 593     aN1V.Subtract(aScaleV * aN0);
 594     aN1V /= aScale;
 595
 596     theD1U += myOffset * aSign * aN1U;
 597     theD1V += myOffset * aSign * aN1V;
 598
 599     return;
 600   }
 601
 602   Standard_Integer OrderU, OrderV;
 603   TColgp_Array2OfVec DerNUV(0, MaxOrder + 1, 0, MaxOrder + 1);
 604   TColgp_Array2OfVec DerSurf(0, MaxOrder + 2, 0, MaxOrder + 2);
 605   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
 606   Bounds(Umin, Umax, Vmin, Vmax);
 607
 608   DerSurf.SetValue(1, 0, theD1U);
 609   DerSurf.SetValue(0, 1, theD1V);
 610   DerSurf.SetValue(1, 1, theD2UV);
 611   DerSurf.SetValue(2, 0, theD2U);
 612   DerSurf.SetValue(0, 2, theD2V);
 613   if (!myBaseSurf.IsNull())
 614     derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
 615   else
 616     derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
 617
 618   gp_Dir Normal;
 619   CSLib_NormalStatus NStatus;
 620   CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax, NStatus, Normal, OrderU, OrderV);
 621   if (NStatus == CSLib_InfinityOfSolutions)
 622   {
 623     gp_Vec aNewDU = theD1U;
 624     gp_Vec aNewDV = theD1V;
 625     // Replace zero derivative and try to calculate normal
 626     if (ReplaceDerivative(theU, theV, aNewDU, aNewDV, the_D1MagTol * the_D1MagTol))
 627     {
 628       DerSurf.SetValue(1, 0, aNewDU);
 629       DerSurf.SetValue(0, 1, aNewDV);
 630       if (!myBaseSurf.IsNull())
 631         derivatives(MaxOrder, 2, theU, theV, myBaseSurf, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
 632       else
 633         derivatives(MaxOrder, 2, theU, theV, myBaseAdaptor, 1, 1, AlongU, AlongV, L, DerNUV, DerSurf);
 634       CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
 635                     NStatus, Normal, OrderU, OrderV);
 636     }
 637   }
 638
 639   if (NStatus != CSLib_Defined)
 640     throw Geom_UndefinedValue(
 641         "GeomEvaluator_OffsetSurface::CalculateD1(): Unable to calculate normal");
 642
 643   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
 644
 645   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
 646   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
 647 }
 648
 649 void GeomEvaluator_OffsetSurface::CalculateD2(
 650     const Standard_Real theU, const Standard_Real theV,
 651     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 652     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
 653     const gp_Vec& theD3U, const gp_Vec& theD3V, const gp_Vec& theD3UUV, const gp_Vec& theD3UVV) const
 654 {
 655   gp_Dir Normal;
 656   CSLib_NormalStatus NStatus;
 657   CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal);
 658
 659   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
 660   Standard_Integer OrderU, OrderV;
 661   TColgp_Array2OfVec DerNUV(0, MaxOrder + 2, 0, MaxOrder + 2);
 662   TColgp_Array2OfVec DerSurf(0, MaxOrder + 3, 0, MaxOrder + 3);
 663
 664   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
 665   Bounds(Umin, Umax, Vmin, Vmax);
 666
 667   DerSurf.SetValue(1, 0, theD1U);
 668   DerSurf.SetValue(0, 1, theD1V);
 669   DerSurf.SetValue(1, 1, theD2UV);
 670   DerSurf.SetValue(2, 0, theD2U);
 671   DerSurf.SetValue(0, 2, theD2V);
 672   DerSurf.SetValue(3, 0, theD3U);
 673   DerSurf.SetValue(2, 1, theD3UUV);
 674   DerSurf.SetValue(1, 2, theD3UVV);
 675   DerSurf.SetValue(0, 3, theD3V);
 676   //*********************
 677
 678   Handle(Geom_BSplineSurface) L;
 679   Standard_Boolean isOpposite = Standard_False;
 680   Standard_Boolean AlongU = Standard_False;
 681   Standard_Boolean AlongV = Standard_False;
 682   if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
 683   {
 684     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
 685     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
 686   }
 687   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
 688
 689   if (!myBaseSurf.IsNull())
 690     derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
 691   else
 692     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 2, 2, AlongU, AlongV, L, DerNUV, DerSurf);
 693
 694   CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
 695                 NStatus, Normal, OrderU, OrderV);
 696   if (NStatus != CSLib_Defined)
 697     throw Geom_UndefinedValue(
 698         "GeomEvaluator_OffsetSurface::CalculateD2(): Unable to calculate normal");
 699
 700   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
 701
 702   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
 703   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
 704
 705   if (!myBaseSurf.IsNull())
 706   {
 707     theD2U = myBaseSurf->DN(theU, theV, 2, 0);
 708     theD2V = myBaseSurf->DN(theU, theV, 0, 2);
 709     theD2UV = myBaseSurf->DN(theU, theV, 1, 1);
 710   }
 711   else
 712   {
 713     theD2U = myBaseAdaptor->DN(theU, theV, 2, 0);
 714     theD2V = myBaseAdaptor->DN(theU, theV, 0, 2);
 715     theD2UV = myBaseAdaptor->DN(theU, theV, 1, 1);
 716   }
 717
 718   theD2U += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV);
 719   theD2V += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV);
 720   theD2UV += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV);
 721 }
 722
 723 void GeomEvaluator_OffsetSurface::CalculateD3(
 724     const Standard_Real theU, const Standard_Real theV,
 725     gp_Pnt& theValue, gp_Vec& theD1U, gp_Vec& theD1V,
 726     gp_Vec& theD2U, gp_Vec& theD2V, gp_Vec& theD2UV,
 727     gp_Vec& theD3U, gp_Vec& theD3V, gp_Vec& theD3UUV, gp_Vec& theD3UVV) const
 728 {
 729   gp_Dir Normal;
 730   CSLib_NormalStatus NStatus;
 731   CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal);
 732   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
 733   Standard_Integer OrderU, OrderV;
 734   TColgp_Array2OfVec DerNUV(0, MaxOrder + 3, 0, MaxOrder + 3);
 735   TColgp_Array2OfVec DerSurf(0, MaxOrder + 4, 0, MaxOrder + 4);
 736   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
 737   Bounds(Umin, Umax, Vmin, Vmax);
 738
 739   DerSurf.SetValue(1, 0, theD1U);
 740   DerSurf.SetValue(0, 1, theD1V);
 741   DerSurf.SetValue(1, 1, theD2UV);
 742   DerSurf.SetValue(2, 0, theD2U);
 743   DerSurf.SetValue(0, 2, theD2V);
 744   DerSurf.SetValue(3, 0, theD3U);
 745   DerSurf.SetValue(2, 1, theD3UUV);
 746   DerSurf.SetValue(1, 2, theD3UVV);
 747   DerSurf.SetValue(0, 3, theD3V);
 748
 749
 750   //*********************
 751   Handle(Geom_BSplineSurface) L;
 752   Standard_Boolean isOpposite = Standard_False;
 753   Standard_Boolean AlongU = Standard_False;
 754   Standard_Boolean AlongV = Standard_False;
 755   if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
 756   {
 757     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
 758     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
 759   }
 760   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
 761
 762   if (!myBaseSurf.IsNull())
 763     derivatives(MaxOrder, 3, theU, theV, myBaseSurf, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
 764   else
 765     derivatives(MaxOrder, 3, theU, theV, myBaseAdaptor, 3, 3, AlongU, AlongV, L, DerNUV, DerSurf);
 766
 767   CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
 768                 NStatus, Normal, OrderU, OrderV);
 769   if (NStatus != CSLib_Defined)
 770     throw Geom_UndefinedValue(
 771         "GeomEvaluator_OffsetSurface::CalculateD3(): Unable to calculate normal");
 772
 773   theValue.SetXYZ(theValue.XYZ() + myOffset * aSign * Normal.XYZ());
 774
 775   theD1U = DerSurf(1, 0) + myOffset * aSign * CSLib::DNNormal(1, 0, DerNUV, OrderU, OrderV);
 776   theD1V = DerSurf(0, 1) + myOffset * aSign * CSLib::DNNormal(0, 1, DerNUV, OrderU, OrderV);
 777
 778   if (!myBaseSurf.IsNull())
 779   {
 780     theD2U   = myBaseSurf->DN(theU, theV, 2, 0);
 781     theD2V   = myBaseSurf->DN(theU, theV, 0, 2);
 782     theD2UV  = myBaseSurf->DN(theU, theV, 1, 1);
 783     theD3U   = myBaseSurf->DN(theU, theV, 3, 0);
 784     theD3V   = myBaseSurf->DN(theU, theV, 0, 3);
 785     theD3UUV = myBaseSurf->DN(theU, theV, 2, 1);
 786     theD3UVV = myBaseSurf->DN(theU, theV, 1, 2);
 787   }
 788   else
 789   {
 790     theD2U   = myBaseAdaptor->DN(theU, theV, 2, 0);
 791     theD2V   = myBaseAdaptor->DN(theU, theV, 0, 2);
 792     theD2UV  = myBaseAdaptor->DN(theU, theV, 1, 1);
 793     theD3U   = myBaseAdaptor->DN(theU, theV, 3, 0);
 794     theD3V   = myBaseAdaptor->DN(theU, theV, 0, 3);
 795     theD3UUV = myBaseAdaptor->DN(theU, theV, 2, 1);
 796     theD3UVV = myBaseAdaptor->DN(theU, theV, 1, 2);
 797   }
 798
 799   theD2U   += aSign * myOffset * CSLib::DNNormal(2, 0, DerNUV, OrderU, OrderV);
 800   theD2V   += aSign * myOffset * CSLib::DNNormal(0, 2, DerNUV, OrderU, OrderV);
 801   theD2UV  += aSign * myOffset * CSLib::DNNormal(1, 1, DerNUV, OrderU, OrderV);
 802   theD3U   += aSign * myOffset * CSLib::DNNormal(3, 0, DerNUV, OrderU, OrderV);
 803   theD3V   += aSign * myOffset * CSLib::DNNormal(0, 3, DerNUV, OrderU, OrderV);
 804   theD3UUV += aSign * myOffset * CSLib::DNNormal(2, 1, DerNUV, OrderU, OrderV);
 805   theD3UVV += aSign * myOffset * CSLib::DNNormal(1, 2, DerNUV, OrderU, OrderV);
 806 }
 807
 808 gp_Vec GeomEvaluator_OffsetSurface::CalculateDN(
 809     const Standard_Real theU, const Standard_Real theV,
 810     const Standard_Integer theNu, const Standard_Integer theNv,
 811     const gp_Vec& theD1U, const gp_Vec& theD1V) const
 812 {
 813   gp_Dir Normal;
 814   CSLib_NormalStatus NStatus;
 815   CSLib::Normal(theD1U, theD1V, the_D1MagTol, NStatus, Normal);
 816   const Standard_Integer MaxOrder = (NStatus == CSLib_Defined) ? 0 : 3;
 817   Standard_Integer OrderU, OrderV;
 818   TColgp_Array2OfVec DerNUV(0, MaxOrder + theNu, 0, MaxOrder + theNv);
 819   TColgp_Array2OfVec DerSurf(0, MaxOrder + theNu + 1, 0, MaxOrder + theNv + 1);
 820
 821   Standard_Real Umin = 0, Umax = 0, Vmin = 0, Vmax = 0;
 822   Bounds(Umin, Umax, Vmin, Vmax);
 823
 824   DerSurf.SetValue(1, 0, theD1U);
 825   DerSurf.SetValue(0, 1, theD1V);
 826
 827   //*********************
 828   Handle(Geom_BSplineSurface) L;
 829   //   Is there any osculatingsurface along U or V;
 830   Standard_Boolean isOpposite = Standard_False;
 831   Standard_Boolean AlongU = Standard_False;
 832   Standard_Boolean AlongV = Standard_False;
 833   if ((NStatus != CSLib_Defined) && !myOscSurf.IsNull())
 834   {
 835     AlongU = myOscSurf->UOscSurf(theU, theV, isOpposite, L);
 836     AlongV = myOscSurf->VOscSurf(theU, theV, isOpposite, L);
 837   }
 838   const Standard_Real aSign = ((AlongV || AlongU) && isOpposite) ? -1. : 1.;
 839
 840   if (!myBaseSurf.IsNull())
 841     derivatives(MaxOrder, 1, theU, theV, myBaseSurf, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
 842   else
 843     derivatives(MaxOrder, 1, theU, theV, myBaseAdaptor, theNu, theNv, AlongU, AlongV, L, DerNUV, DerSurf);
 844
 845   CSLib::Normal(MaxOrder, DerNUV, the_D1MagTol, theU, theV, Umin, Umax, Vmin, Vmax,
 846                 NStatus, Normal, OrderU, OrderV);
 847   if (NStatus != CSLib_Defined)
 848     throw Geom_UndefinedValue(
 849         "GeomEvaluator_OffsetSurface::CalculateDN(): Unable to calculate normal");
 850
 851   gp_Vec D;
 852   if (!myBaseSurf.IsNull())
 853     D = myBaseSurf->DN(theU, theV, theNu, theNv);
 854   else
 855     D = myBaseAdaptor->DN(theU, theV, theNu, theNv);
 856
 857   D += aSign * myOffset * CSLib::DNNormal(theNu, theNv, DerNUV, OrderU, OrderV);
 858   return D;
 859 }
 860
 861
 862 void GeomEvaluator_OffsetSurface::Bounds(Standard_Real& theUMin, Standard_Real& theUMax,
 863                                          Standard_Real& theVMin, Standard_Real& theVMax) const
 864 {
 865   if (!myBaseSurf.IsNull())
 866     myBaseSurf->Bounds(theUMin, theUMax, theVMin, theVMax);
 867   else
 868   {
 869     theUMin = myBaseAdaptor->FirstUParameter();
 870     theUMax = myBaseAdaptor->LastUParameter();
 871     theVMin = myBaseAdaptor->FirstVParameter();
 872     theVMax = myBaseAdaptor->LastVParameter();
 873   }
 874 }
 875
 876
 877 Standard_Boolean GeomEvaluator_OffsetSurface::ReplaceDerivative(
 878     const Standard_Real theU, const Standard_Real theV,
 879     gp_Vec& theDU, gp_Vec& theDV,
 880     const Standard_Real theSquareTol) const
 881 {
 882   Standard_Boolean isReplaceDU = theDU.SquareMagnitude() < theSquareTol;
 883   Standard_Boolean isReplaceDV = theDV.SquareMagnitude() < theSquareTol;
 884   Standard_Boolean isReplaced = Standard_False;
 885   if (isReplaceDU ^ isReplaceDV)
 886   {
 887     Standard_Real aU = theU;
 888     Standard_Real aV = theV;
 889     Standard_Real aUMin = 0, aUMax = 0, aVMin = 0, aVMax = 0;
 890     Bounds(aUMin, aUMax, aVMin, aVMax);
 891
 892     // Calculate step along non-zero derivative
 893     Standard_Real aStep;
 894     Handle(Adaptor3d_HSurface) aSurfAdapt;
 895     if (!myBaseAdaptor.IsNull())
 896       aSurfAdapt = myBaseAdaptor;
 897     else
 898       aSurfAdapt = new GeomAdaptor_HSurface(myBaseSurf);
 899     if (isReplaceDV)
 900     {
 901       aStep = Precision::Confusion() * theDU.Magnitude();
 902       if (aStep > aUMax - aUMin)
 903         aStep = (aUMax - aUMin) / 100.;
 904     }
 905     else
 906     {
 907       aStep = Precision::Confusion() * theDV.Magnitude();
 908       if (aStep > aVMax - aVMin)
 909         aStep = (aVMax - aVMin) / 100.;
 910     }
 911
 912     gp_Pnt aP;
 913     gp_Vec aDU, aDV;
 914     // Step away from currect parametric coordinates and calculate derivatives once again.
 915     // Replace zero derivative by the obtained.
 916     for (Standard_Real aStepSign = -1.0; aStepSign <= 1.0 && !isReplaced; aStepSign += 2.0)
 917     {
 918       if (isReplaceDV)
 919       {
 920         aU = theU + aStepSign * aStep;
 921         if (aU < aUMin || aU > aUMax)
 922           continue;
 923       }
 924       else
 925       {
 926         aV = theV + aStepSign * aStep;
 927         if (aV < aVMin || aV > aVMax)
 928           continue;
 929       }
 930
 931       BaseD1(aU, aV, aP, aDU, aDV);
 932
 933       if (isReplaceDU && aDU.SquareMagnitude() > theSquareTol)
 934       {
 935         theDU = aDU;
 936         isReplaced = Standard_True;
 937       }
 938       if (isReplaceDV && aDV.SquareMagnitude() > theSquareTol)
 939       {
 940         theDV = aDV;
 941         isReplaced = Standard_True;
 942       }
 943     }
 944   }
 945   return isReplaced;
 946 }