0028838: Configuration - undefine macros coming from X11 headers in place of collision
[occt.git] / src / LProp / LProp_SLProps.gxx
1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 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 <LProp_Status.hxx>
16 #include <LProp_NotDefined.hxx>
17 #include <Standard_OutOfRange.hxx>
18 #include <Standard_DomainError.hxx>
19 #include <CSLib.hxx>
20 #include <CSLib_DerivativeStatus.hxx>
21 #include <CSLib_NormalStatus.hxx>
22 #include <TColgp_Array2OfVec.hxx>
23 #include <math_DirectPolynomialRoots.hxx>
24
25 static const Standard_Real MinStep   = 1.0e-7;
26
27 static Standard_Boolean IsTangentDefined (LProp_SLProps& SProp,
28                                           const Standard_Integer cn,
29                                           const Standard_Real linTol,
30                                           const Standard_Integer  Derivative, 
31                                           Standard_Integer&       Order,
32                                           LProp_Status&  theStatus)
33 {
34   Standard_Real Tol = linTol * linTol;
35   gp_Vec V[2];
36   Order = 0;
37
38   while (Order < 3)
39   {
40     Order++;
41     if(cn >= Order)
42     {
43       switch(Order)
44       {
45       case 1:
46         V[0] = SProp.D1U();
47         V[1] = SProp.D1V();
48         break;
49       case 2:
50         V[0] = SProp.D2U();
51         V[1] = SProp.D2V();
52         break;
53       }//switch(Order)
54
55       if(V[Derivative].SquareMagnitude() > Tol)
56       {
57         theStatus = LProp_Defined;
58         return Standard_True;
59       }
60     }//if(cn >= Order)
61     else
62     {
63       theStatus = LProp_Undefined;
64       return Standard_False;
65     }
66   }
67
68   return Standard_False;
69 }
70
71 LProp_SLProps::LProp_SLProps (const Surface& S, 
72                               const Standard_Real U,  
73                               const Standard_Real V, 
74                               const Standard_Integer N, 
75                               const Standard_Real Resolution) 
76         : mySurf(S),myDerOrder(N), myCN(4), // (Tool::Continuity(S)),
77             myLinTol(Resolution)
78 {
79   Standard_OutOfRange_Raise_if(N < 0 || N > 2, 
80                         "LProp_SLProps::LProp_SLProps()");
81
82   SetParameters(U, V);
83 }
84
85 LProp_SLProps::LProp_SLProps (const Surface& S, 
86                               const Standard_Integer  N, 
87                               const Standard_Real     Resolution) 
88         : mySurf(S), myU(RealLast()), myV(RealLast()), myDerOrder(N),
89             myCN(4), // (Tool::Continuity(S))
90             myLinTol(Resolution), 
91             myUTangentStatus (LProp_Undecided),
92             myVTangentStatus (LProp_Undecided),
93             myNormalStatus   (LProp_Undecided),
94             myCurvatureStatus(LProp_Undecided)
95 {
96   Standard_OutOfRange_Raise_if(N < 0 || N > 2, 
97                       "LProp_SLProps::LProp_SLProps()");
98 }
99
100 LProp_SLProps::LProp_SLProps (const Standard_Integer  N, 
101                               const Standard_Real Resolution) 
102         : myU(RealLast()), myV(RealLast()), myDerOrder(N), myCN(0),
103         myLinTol(Resolution),
104         myUTangentStatus (LProp_Undecided),
105         myVTangentStatus (LProp_Undecided),
106         myNormalStatus   (LProp_Undecided),
107         myCurvatureStatus(LProp_Undecided)
108 {
109   Standard_OutOfRange_Raise_if(N < 0 || N > 2, 
110                 "LProp_SLProps::LProp_SLProps() bad level");
111 }
112
113 void LProp_SLProps::SetSurface (const Surface& S ) 
114 {
115   mySurf = S;
116   myCN = 4; // =Tool::Continuity(S);
117 }
118
119 void LProp_SLProps::SetParameters (const Standard_Real U, 
120                                    const Standard_Real V)
121 {
122   myU = U;
123   myV = V;
124   switch (myDerOrder)
125   {
126   case 0:
127     Tool::Value(mySurf, myU, myV, myPnt);
128     break;
129   case 1:
130     Tool::D1(mySurf, myU, myV, myPnt, myD1u, myD1v);
131     break;
132   case 2:
133     Tool::D2(mySurf, myU, myV, myPnt, myD1u, myD1v, myD2u, myD2v, myDuv);
134     break;
135   }
136
137   myUTangentStatus  = LProp_Undecided;
138   myVTangentStatus  = LProp_Undecided;
139   myNormalStatus    = LProp_Undecided;
140   myCurvatureStatus = LProp_Undecided;
141 }
142
143 const gp_Pnt& LProp_SLProps::Value() const
144 {
145   return myPnt;
146 }
147
148 const gp_Vec& LProp_SLProps::D1U()
149 {
150   if (myDerOrder < 1)
151   {
152     myDerOrder =1;
153     Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v);
154   }
155
156   return myD1u;
157 }
158
159 const gp_Vec& LProp_SLProps::D1V()
160 {
161   if (myDerOrder < 1)
162   {
163     myDerOrder =1;
164     Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v);
165   }
166
167   return myD1v;
168 }
169
170 const gp_Vec& LProp_SLProps::D2U()
171 {
172   if (myDerOrder < 2) 
173   {
174     myDerOrder =2;
175     Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
176   }
177
178   return myD2u;
179 }
180
181 const gp_Vec& LProp_SLProps::D2V()
182 {
183   if (myDerOrder < 2) 
184   {
185     myDerOrder =2;
186     Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
187   }
188
189   return myD2v;
190 }
191
192 const gp_Vec& LProp_SLProps::DUV()
193 {
194   if (myDerOrder < 2) 
195   {
196     myDerOrder =2;
197     Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
198   }
199
200   return myDuv;
201 }
202
203 Standard_Boolean LProp_SLProps::IsTangentUDefined ()
204 {
205   if (myUTangentStatus == LProp_Undefined)
206     return Standard_False;
207   else if (myUTangentStatus >= LProp_Defined)
208     return Standard_True; 
209
210   // uTangentStatus == Lprop_Undecided 
211   // we have to calculate the first non null U derivative
212   return IsTangentDefined(*this, myCN, myLinTol, 0, 
213           mySignificantFirstDerivativeOrderU, myUTangentStatus);
214 }
215
216 void LProp_SLProps::TangentU (gp_Dir& D) 
217 {
218   if(!IsTangentUDefined())
219     throw LProp_NotDefined();
220
221   if(mySignificantFirstDerivativeOrderU == 1)
222     D = gp_Dir(myD1u);
223   else
224   {
225     const Standard_Real DivisionFactor = 1.e-3;
226     Standard_Real anUsupremum, anUinfium;
227     Standard_Real anVsupremum, anVinfium;
228     Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum);
229     Standard_Real du;
230     if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
231       du = 0.0;
232     else
233       du = anUsupremum-anUinfium;
234     
235     const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);
236
237     gp_Vec V = myD2u;
238     
239     Standard_Real u;
240           
241     if(myU-anUinfium < aDeltaU)
242       u = myU+aDeltaU;
243     else
244       u = myU-aDeltaU;
245     
246     gp_Pnt P1, P2;
247     Tool::Value(mySurf, Min(myU, u),myV,P1);
248     Tool::Value(mySurf, Max(myU, u),myV,P2);
249     
250     gp_Vec V1(P1,P2);
251     Standard_Real aDirFactor = V.Dot(V1);
252     
253     if(aDirFactor < 0.0)
254       V = -V;
255     
256     D = gp_Dir(V);
257   }
258 }
259
260 Standard_Boolean LProp_SLProps::IsTangentVDefined ()
261 {
262   if (myVTangentStatus == LProp_Undefined)
263     return Standard_False;
264   else if (myVTangentStatus >= LProp_Defined)
265     return Standard_True; 
266
267   // vTangentStatus == Lprop_Undecided 
268   // we have to calculate the first non null V derivative
269   return IsTangentDefined(*this, myCN, myLinTol, 1, 
270           mySignificantFirstDerivativeOrderV, myVTangentStatus);
271 }
272
273 void LProp_SLProps::TangentV (gp_Dir& D)
274 {
275   if(!IsTangentVDefined())
276     throw LProp_NotDefined();
277   
278   if(mySignificantFirstDerivativeOrderV == 1)
279     D = gp_Dir(myD1v);
280   else
281   {
282     const Standard_Real DivisionFactor = 1.e-3;
283     Standard_Real anUsupremum, anUinfium;
284     Standard_Real anVsupremum, anVinfium;
285     Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum);
286     
287     Standard_Real dv;
288     if((anVsupremum >= RealLast()) || (anVinfium <= RealFirst()))
289       dv = 0.0;
290     else
291       dv = anVsupremum-anVinfium;
292     
293     const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep);
294
295     gp_Vec V = myD2v;
296     
297     Standard_Real v;
298           
299     if(myV-anVinfium < aDeltaV)
300       v = myV+aDeltaV;
301     else
302       v = myV-aDeltaV;
303     
304     gp_Pnt P1, P2;
305     Tool::Value(mySurf, myU, Min(myV, v),P1);
306     Tool::Value(mySurf, myU, Max(myV, v),P2);
307     
308     gp_Vec V1(P1,P2);
309     Standard_Real aDirFactor = V.Dot(V1);
310     
311     if(aDirFactor < 0.0)
312       V = -V;
313     
314     D = gp_Dir(V);
315   }
316 }
317
318 Standard_Boolean LProp_SLProps::IsNormalDefined()
319 {
320   if (myNormalStatus == LProp_Undefined)
321     return Standard_False;
322   else if (myNormalStatus >= LProp_Defined)
323     return Standard_True;
324
325   // status = UnDecided 
326   // first try the standard computation of the normal.
327   CSLib_DerivativeStatus aStatus = CSLib_Done;
328   CSLib::Normal(myD1u, myD1v, myLinTol, aStatus, myNormal);
329   if (aStatus == CSLib_Done)
330   {
331     myNormalStatus = LProp_Computed;
332     return Standard_True;
333   }
334
335   // else solve the degenerated case only if continuity >= 2
336
337   myNormalStatus = LProp_Undefined;
338   return Standard_False;
339 }
340
341 const gp_Dir& LProp_SLProps::Normal ()
342 {
343   if(!IsNormalDefined())
344   {
345     throw LProp_NotDefined();
346   }
347
348   return myNormal;
349 }
350
351
352 Standard_Boolean LProp_SLProps::IsCurvatureDefined ()
353 {
354   if (myCurvatureStatus == LProp_Undefined)
355     return Standard_False;
356   else if (myCurvatureStatus >= LProp_Defined)
357     return Standard_True;
358   
359   if(myCN < 2)
360   {
361     myCurvatureStatus = LProp_Undefined;
362     return Standard_False;
363   }
364
365   // status = UnDecided 
366   if (!IsNormalDefined())
367   {
368     myCurvatureStatus = LProp_Undefined;
369     return Standard_False;
370   }
371
372   // pour eviter un plantage dans le cas du caro pointu
373   // en fait on doit pouvoir calculer les courbure
374   // avoir
375   if(!IsTangentUDefined() || !IsTangentVDefined())
376   {
377     myCurvatureStatus = LProp_Undefined;
378     return Standard_False;
379   }
380
381   // here we compute the curvature features of the surface
382
383   gp_Vec Norm (myNormal);
384
385   Standard_Real E = myD1u.SquareMagnitude();
386   Standard_Real F = myD1u.Dot(myD1v);
387   Standard_Real G = myD1v.SquareMagnitude();
388
389   if(myDerOrder < 2)
390     this->D2U();
391
392   Standard_Real L = Norm.Dot(myD2u);
393   Standard_Real M = Norm.Dot(myDuv);
394   Standard_Real N = Norm.Dot(myD2v);
395
396   Standard_Real A = E * M - F * L;
397   Standard_Real B = E * N - G * L;
398   Standard_Real C = F * N - G * M;
399
400   Standard_Real MaxABC = Max(Max(Abs(A),Abs(B)),Abs(C));
401   if (MaxABC < RealEpsilon())    // ombilic
402   {
403     myMinCurv = N / G;
404     myMaxCurv = myMinCurv;
405     myDirMinCurv = gp_Dir (myD1u);
406     myDirMaxCurv = gp_Dir (myD1u.Crossed(Norm));
407     myMeanCurv = myMinCurv;            // (Cmin + Cmax) / 2.
408     myGausCurv = myMinCurv * myMinCurv;  // (Cmin * Cmax)
409     myCurvatureStatus = LProp_Computed;
410     return Standard_True;
411   }
412
413   A = A / MaxABC;
414   B = B / MaxABC;
415   C = C / MaxABC;
416   Standard_Real Curv1, Curv2, Root1, Root2;
417   gp_Vec VectCurv1, VectCurv2;
418
419   if (Abs(A) > RealEpsilon())
420   {
421     math_DirectPolynomialRoots Root (A, B, C);
422     if(Root.NbSolutions() != 2)
423     {
424       myCurvatureStatus = LProp_Undefined;
425       return Standard_False;
426     }
427     else
428     {
429       Root1 = Root.Value(1);
430       Root2 = Root.Value(2);
431       Curv1 = ((L * Root1 + 2. * M) * Root1 + N) /
432               ((E * Root1 + 2. * F) * Root1 + G);
433       Curv2 = ((L * Root2 + 2. * M) * Root2 + N) /
434               ((E * Root2 + 2. * F) * Root2 + G);
435       VectCurv1 = Root1 * myD1u + myD1v;
436       VectCurv2 = Root2 * myD1u + myD1v;
437     }
438   }
439   else if (Abs(C) > RealEpsilon())
440   {
441     math_DirectPolynomialRoots Root(C, B, A);
442
443     if((Root.NbSolutions() != 2))
444     {
445       myCurvatureStatus = LProp_Undefined;
446       return Standard_False;
447     }
448     else
449     {
450       Root1 = Root.Value(1);
451       Root2 = Root.Value(2);
452       Curv1 = ((N * Root1 + 2. * M) * Root1 + L) /
453               ((G * Root1 + 2. * F) * Root1 + E);
454       Curv2 = ((N * Root2 + 2. * M) * Root2 + L) /
455               ((G * Root2 + 2. * F) * Root2 + E);
456       VectCurv1 = myD1u + Root1 * myD1v;
457       VectCurv2 = myD1u + Root2 * myD1v;
458     }
459   }
460   else
461   {
462     Curv1 = L / E;
463     Curv2 = N / G;
464     VectCurv1 = myD1u;
465     VectCurv2 = myD1v;
466   }
467
468   if (Curv1 < Curv2)
469   {
470     myMinCurv = Curv1;
471     myMaxCurv = Curv2;
472     myDirMinCurv = gp_Dir (VectCurv1);
473     myDirMaxCurv = gp_Dir (VectCurv2);
474   }
475   else
476   {
477     myMinCurv = Curv2;
478     myMaxCurv = Curv1;
479     myDirMinCurv = gp_Dir (VectCurv2);
480     myDirMaxCurv = gp_Dir (VectCurv1);
481   }
482
483   myMeanCurv = ((N * E) - (2. * M * F) + (L * G)) // voir Farin p.282
484                / (2. * ((E * G) - (F * F)));
485   myGausCurv = ((L * N) - (M * M)) 
486                / ((E * G) - (F * F));
487   myCurvatureStatus = LProp_Computed;
488   return Standard_True;
489 }
490
491 Standard_Boolean LProp_SLProps::IsUmbilic ()
492 {
493   if(!IsCurvatureDefined())
494     throw LProp_NotDefined();
495   
496   return Abs(myMaxCurv - myMinCurv) < Abs(Epsilon(myMaxCurv));
497 }
498
499 Standard_Real LProp_SLProps::MaxCurvature ()
500 {
501   if(!IsCurvatureDefined())
502     throw LProp_NotDefined();
503
504   return myMaxCurv;
505 }
506
507 Standard_Real LProp_SLProps::MinCurvature ()
508 {
509   if(!IsCurvatureDefined())
510     throw LProp_NotDefined();
511   
512   return myMinCurv;
513 }
514
515 void LProp_SLProps::CurvatureDirections(gp_Dir& Max, gp_Dir& Min)
516 {
517   if(!IsCurvatureDefined())
518     throw LProp_NotDefined();
519
520   Max = myDirMaxCurv;
521   Min = myDirMinCurv;
522 }
523
524 Standard_Real LProp_SLProps::MeanCurvature ()
525 {
526   if(!IsCurvatureDefined())
527     throw LProp_NotDefined();
528
529   return myMeanCurv;
530 }
531
532 Standard_Real LProp_SLProps::GaussianCurvature ()
533 {
534   if(!IsCurvatureDefined())
535     throw LProp_NotDefined();
536
537   return myGausCurv;
538 }