d12bd588c706d6d9832ed37de0b893d89554913a
[occt.git] / 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-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21
22 #include <GCPnts_DeflectionType.hxx>
23 #include <Standard_ConstructionError.hxx>
24 #include <Precision.hxx>
25 #include <gp_XYZ.hxx>
26 #include <gp_Pnt.hxx>
27 #include <gp_Vec.hxx>
28 #include <gp.hxx>
29
30 #define Us3 0.3333333333333333333333333333
31
32 void GCPnts_TangentialDeflection::EvaluateDu (
33   const TheCurve& C,
34   const Standard_Real      U,
35   gp_Pnt&   P,
36   Standard_Real&     Du,
37   Standard_Boolean&  NotDone) const {
38
39   gp_Vec T, N;
40   D2 (C, U, P, T, N);
41   Standard_Real Lt   = T.Magnitude ();
42   Standard_Real LTol = Precision::Confusion ();
43   if (Lt > LTol && N.Magnitude () > LTol) {
44     Standard_Real Lc = N.CrossMagnitude (T);
45     Standard_Real Ln = Lc/Lt;
46     if (Ln > LTol) {
47       Du = sqrt (8.0 * curvatureDeflection / Ln);
48       NotDone = Standard_False;
49     }
50   }
51 }
52
53
54 //=======================================================================
55 //function : GCPnts_TangentialDeflection
56 //purpose  : 
57 //=======================================================================
58
59 GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
60  const TheCurve&        C,
61  const Standard_Real    AngularDeflection,
62  const Standard_Real    CurvatureDeflection,
63  const Standard_Integer MinimumOfPoints,
64  const Standard_Real    UTol)
65
66
67   Initialize (C,AngularDeflection,CurvatureDeflection,MinimumOfPoints,UTol); 
68 }
69
70
71 //=======================================================================
72 //function : GCPnts_TangentialDeflection
73 //purpose  : 
74 //=======================================================================
75
76 GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
77  const TheCurve&        C,
78  const Standard_Real    FirstParameter,
79  const Standard_Real    LastParameter,
80  const Standard_Real    AngularDeflection,
81  const Standard_Real    CurvatureDeflection,
82  const Standard_Integer MinimumOfPoints,
83  const Standard_Real    UTol)
84
85
86   Initialize (C, 
87         FirstParameter, 
88         LastParameter,
89         AngularDeflection, 
90         CurvatureDeflection, 
91         MinimumOfPoints,
92         UTol);
93 }
94
95
96
97 //=======================================================================
98 //function : Initialize
99 //purpose  : 
100 //=======================================================================
101
102 void GCPnts_TangentialDeflection::Initialize (
103  const TheCurve&        C, 
104  const Standard_Real    AngularDeflection, 
105  const Standard_Real    CurvatureDeflection,
106  const Standard_Integer MinimumOfPoints,
107  const Standard_Real    UTol)
108
109 {
110   Initialize (C, 
111               C.FirstParameter (), 
112               C.LastParameter (),
113               AngularDeflection, 
114               CurvatureDeflection,
115               MinimumOfPoints,
116               UTol);
117 }
118
119
120 //=======================================================================
121 //function : Initialize
122 //purpose  : 
123 //=======================================================================
124
125 void GCPnts_TangentialDeflection::Initialize (
126  const TheCurve&        C, 
127  const Standard_Real    FirstParameter, 
128  const Standard_Real    LastParameter,
129  const Standard_Real    AngularDeflection, 
130  const Standard_Real    CurvatureDeflection,
131  const Standard_Integer MinimumOfPoints,
132  const Standard_Real    UTol)
133
134 {
135   
136   Standard_ConstructionError_Raise_if (CurvatureDeflection <= Precision::Confusion () || AngularDeflection   <= Precision::Angular (), "GCPnts_TangentialDeflection::Initialize - Zero Deflection")
137
138  parameters.Clear ();
139  points    .Clear ();
140  if (FirstParameter < LastParameter) {
141    firstu = FirstParameter;
142    lastu  = LastParameter;
143  }
144  else {
145    lastu  = FirstParameter;
146    firstu = LastParameter;
147  }
148  uTol                = UTol;
149  angularDeflection   = AngularDeflection;
150  curvatureDeflection = CurvatureDeflection;
151  minNbPnts           = Max (MinimumOfPoints, 2);
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 //=======================================================================
211 //function : PerformCircular
212 //purpose  : 
213 //=======================================================================
214
215 void GCPnts_TangentialDeflection::PerformCircular (const TheCurve& C) 
216 {
217   // akm 8/01/02 : check the radius before divide by it
218   Standard_Real dfR = C.Circle().Radius();
219   Standard_Real Du = 0.;
220   if (Abs(dfR) > Precision::Confusion())
221     Du = Max(1.0e0 - (curvatureDeflection/dfR),0.0e0) ;
222   Du  = acos (Du); Du+=Du;
223   Du               = Min (Du, angularDeflection);
224   Standard_Integer NbPoints = (Standard_Integer )((lastu - firstu) / Du);
225   NbPoints         = Max (NbPoints, minNbPnts-1);
226   Du               = (lastu - firstu) / NbPoints;
227
228   gp_Pnt P;
229   Standard_Real U = firstu;
230   for (Standard_Integer i = 1; i <= NbPoints; i++) {
231     D0 (C, U, P);
232     parameters.Append (U);
233     points    .Append (P);
234     U += Du;
235   }
236   D0 (C, lastu, P);
237   parameters.Append (lastu);
238   points    .Append (P);
239 }
240
241
242
243 //=======================================================================
244 //function : PerformCurve
245 //purpose  : On respecte ll'angle et la fleche, on peut imposer un nombre
246 //           minimum de points sur un element lineaire
247 //=======================================================================
248
249 void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
250
251 {
252   Standard_Integer i;       
253   gp_XYZ V1, V2;
254   gp_Pnt MiddlePoint, CurrentPoint, LastPoint;   
255   Standard_Real Du, Dusave, MiddleU, L1, L2;
256
257   Standard_Real U1       = firstu;   
258   Standard_Real LTol     = Precision::Confusion ();  //protection longueur nulle
259   Standard_Real ATol     = Precision::Angular ();    //protection angle nul
260
261   D0 (C, lastu, LastPoint);
262
263   //Initialization du calcul
264
265   Standard_Boolean NotDone = Standard_True;
266   Dusave = (lastu - firstu)*Us3;
267   Du     = Dusave;
268   EvaluateDu (C, U1, CurrentPoint, Du, NotDone);
269   parameters.Append (U1);
270   points    .Append (CurrentPoint);
271
272   if (NotDone) {
273     //C'est soit une droite, soit une singularite :
274     V1 = LastPoint.XYZ ();
275     V1.Subtract (CurrentPoint.XYZ());
276     L1 = V1.Modulus ();
277     if (L1 > LTol) {
278       //Si c'est une droite on verifie en calculant minNbPoints :
279       Standard_Boolean IsLine   = Standard_True;
280       Standard_Integer NbPoints = 3;
281       if (minNbPnts > 3) NbPoints = minNbPnts;
282       Du = (lastu-firstu)/NbPoints;
283       MiddleU = firstu + Du;
284       for (i = 2; i < NbPoints; i++) {
285         D0 (C, MiddleU, MiddlePoint);
286         V2 = MiddlePoint.XYZ();
287         V2.Subtract (CurrentPoint.XYZ());
288         L2 = V2.Modulus ();
289         if (L2 > LTol) {
290           if (((V2.CrossMagnitude (V1))/(L1*L2)) >= ATol) {
291             //C'etait une singularite
292             IsLine = Standard_False;
293             break;
294           }
295           if (minNbPnts > 2) {
296             parameters.Append (MiddleU);
297             points    .Append (MiddlePoint);
298           }
299         }
300         MiddleU += Du;
301       }
302       if (IsLine) {
303         //C'etait une droite (plusieurs poles alignes), Calcul termine :
304         parameters.Append (lastu);
305         points    .Append (LastPoint);
306         return;
307       }
308       else {
309         //c'etait une singularite on continue :
310         Standard_Integer pointsLength=points.Length ();
311         for (i = 2; i <= pointsLength; i++) {
312           points    .Remove (i);
313           parameters.Remove (i);
314           pointsLength--;
315         }
316         Du = Dusave;
317       }
318     }
319     else  {
320       
321       Du = (lastu-firstu)/2.1;
322       MiddleU = firstu + Du;
323       D0 (C, MiddleU, MiddlePoint);
324       V1 = MiddlePoint.XYZ ();
325       V1.Subtract (CurrentPoint.XYZ());
326       L1 = V1.Modulus ();
327       if (L1 < LTol) {
328         // L1 < LTol C'est une courbe de longueur nulle, calcul termine :
329         // on renvoi un segment de 2 points   (protection)
330         parameters.Append (lastu);
331         points    .Append (LastPoint);
332         return;
333       }
334     }
335   }
336
337   if (Du > Dusave) Du = Dusave;
338   else             Dusave = Du;
339
340   if (Du < uTol) {
341     Du = lastu - firstu;
342     if (Du < uTol) {
343       parameters.Append (lastu);
344       points    .Append (LastPoint);
345       return;
346     }
347   }
348
349   //Traitement normal pour une courbe
350   Standard_Boolean MorePoints = Standard_True;
351   Standard_Real U2            = firstu;   
352   Standard_Real AngleMax      = angularDeflection * 0.5;  //car on prend le point milieu
353
354   gp_Pnt aPrevPoint = points.Last();
355
356   while (MorePoints) {
357
358     U2 += Du;                             
359
360     if (U2 >= lastu) {                       //Bout de courbe
361       U2 = lastu;
362       CurrentPoint = LastPoint;
363       Du = U2-U1;
364       Dusave = Du;
365     }
366     else D0 (C, U2, CurrentPoint);           //Point suivant
367
368     Standard_Real Coef, ACoef, FCoef;
369     Standard_Boolean Correction, TooLarge, TooSmall;
370     TooLarge   = Standard_False;
371     TooSmall   = Standard_False;
372     Correction = Standard_True;
373
374     Standard_Real lastCoef = 0;
375      
376     while (Correction) {                     //Ajustement Du
377       MiddleU = (U1+U2)*0.5;                 //Verif / au point milieu
378       D0 (C, MiddleU, MiddlePoint);
379
380       V1 = CurrentPoint.XYZ ();              //Critere de fleche
381       V1.Subtract (aPrevPoint.XYZ());
382       V2 = MiddlePoint.XYZ ();
383       V2.Subtract (aPrevPoint.XYZ());
384       L1 = V1.Modulus ();
385       if (L1 > LTol) FCoef = V1.CrossMagnitude(V2)/(L1*curvatureDeflection);
386       else           FCoef = 0.0;
387
388       V1 = CurrentPoint.XYZ ();              //Critere d'angle
389       V1.Subtract (MiddlePoint.XYZ ());
390       L1 = V1.Modulus ();
391       L2 = V2.Modulus ();
392       Standard_Real angg = V1.CrossMagnitude(V2)/(L1*L2);
393       if (L1 > LTol && L2 > LTol) ACoef = angg/AngleMax;
394       else                        ACoef = 0.0;
395
396       if (ACoef >= FCoef) Coef = ACoef;      //On retient le plus penalisant
397       else                Coef = FCoef;
398
399
400       if (Coef <= 1.0) {
401         if (Abs (lastu-U2) < uTol) {
402           parameters.Append (lastu);
403           points    .Append (LastPoint);
404           MorePoints = Standard_False;
405           Correction = Standard_False;
406         }
407         else {
408           if (Coef >= 0.55 || TooLarge) { 
409             parameters.Append (U2);
410             points    .Append (CurrentPoint);
411             aPrevPoint = CurrentPoint;
412             Correction = Standard_False;
413           }
414           else if (TooSmall) {
415             Correction = Standard_False;
416             aPrevPoint = CurrentPoint;
417           }
418           else {
419             TooSmall = Standard_True;
420             lastCoef = Coef;
421             //Standard_Real UUU2 = U2;
422             Du += Min((U2-U1)*(1.-Coef), Du*Us3);
423
424             U2 = U1 + Du;
425             //if (U2 >= lastu) U2 = UUU2;
426             if (U2 >= lastu) {
427               parameters.Append (lastu);
428               points    .Append (LastPoint);
429               MorePoints = Standard_False;
430               Correction = Standard_False;
431             }
432             else D0 (C, U2, CurrentPoint);
433           }
434         }
435       }
436       else {
437
438         if (Coef >= 1.5) {
439           U2 = MiddleU;
440           Du  = U2-U1;
441           CurrentPoint = MiddlePoint;
442         }
443         else {
444           Du*=0.9;
445           U2 = U1 + Du;
446           D0 (C, U2, CurrentPoint);
447           TooLarge = Standard_True;
448         }
449
450       }
451     }
452     
453     Du  = U2-U1;
454
455     if (MorePoints) {
456       if (U1 > firstu) {
457         if (FCoef > ACoef) {
458           //La fleche est critere de decoupage
459           EvaluateDu (C, U2, CurrentPoint, Du, NotDone);
460           if (NotDone) {
461             Du += (Du-Dusave)*(Du/Dusave);
462             if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
463             if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
464           }
465         }
466         else {
467           //L'angle est le critere de decoupage
468           Du += (Du-Dusave)*(Du/Dusave);
469           if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
470           if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
471         }
472       }
473       
474       if (Du < uTol) {
475         Du = lastu - U2;
476         if (Du < uTol) {
477           parameters.Append (lastu);
478           points    .Append (LastPoint);
479           MorePoints = Standard_False;
480         }
481         else if (Du*Us3 > uTol) Du*=Us3;
482       }
483       U1 = U2;
484       Dusave = Du;
485     }
486   }
487   //Recalage avant dernier point :
488   i = points.Length()-1;
489 //  Real d = points (i).Distance (points (i+1));
490 // if (Abs(parameters (i) - parameters (i+1))<= 0.000001 || d < Precision::Confusion()) {
491 //    cout<<"deux points confondus"<<endl;
492 //    parameters.Remove (i+1);
493 //    points.Remove (i+1);
494 //    i--;
495 //  }
496   
497   if (i >= 2) {
498     MiddleU = parameters (i-1);
499     MiddleU = (lastu + MiddleU)*0.5;
500     D0 (C, MiddleU, MiddlePoint);
501     parameters.SetValue (i, MiddleU);
502     points    .SetValue (i, MiddlePoint);
503   }
504
505   //-- On rajoute des points aux milieux des segments si le nombre
506   //-- mini de points n'est pas atteint
507   //--
508   Standard_Integer Nbp =  points.Length();
509   Standard_Integer MinNb= (9*minNbPnts)/10;
510   //if(MinNb<4) MinNb=4;
511
512   //-- if(Nbp <  MinNb) { cout<<"\n*"; } else {  cout<<"\n."; } 
513   while(Nbp < MinNb) { 
514     //-- cout<<" \nGCPnts TangentialDeflection : Ajout de Points ("<<Nbp<<" "<<minNbPnts<<" )"<<endl;
515     for(i=2; i<=Nbp; i++) { 
516       MiddleU = (parameters.Value(i-1)+parameters.Value(i))*0.5;
517       D0 (C, MiddleU, MiddlePoint); 
518       parameters.InsertBefore(i,MiddleU);
519       points.InsertBefore(i,MiddlePoint);
520       Nbp++;
521       i++;
522     }
523   }
524 }