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