74afe032166e373345f90d0c4b346d966e664bc7
[occt.git] / src / GCPnts / GCPnts_TangentialDeflection.gxx
1 // File:        GCPnts_TangentialDeflection.gxx
2 // Created:     Fri Nov  8 11:26:11 1996
3 // Author:      Jean Claude VAUTHIER
4 //              <jcv@brunox.paris1.matra-dtv.fr>
5
6 #include <GCPnts_DeflectionType.hxx>
7 #include <Standard_ConstructionError.hxx>
8 #include <Precision.hxx>
9 #include <gp_XYZ.hxx>
10 #include <gp_Pnt.hxx>
11 #include <gp_Vec.hxx>
12 #include <gp.hxx>
13
14 #define Us3 0.3333333333333333333333333333
15
16 void GCPnts_TangentialDeflection::EvaluateDu (
17   const TheCurve& C,
18   const Standard_Real      U,
19         gp_Pnt&   P,
20         Standard_Real&     Du,
21         Standard_Boolean&  NotDone) const {
22
23   gp_Vec T, N;
24   D2 (C, U, P, T, N);
25   Standard_Real Lt   = T.Magnitude ();
26   Standard_Real LTol = Precision::Confusion ();
27   if (Lt > LTol && N.Magnitude () > LTol) {
28     Standard_Real Lc = N.CrossMagnitude (T);
29     if (Lc/Lt > LTol) {
30       Du = sqrt (8.0 * curvatureDeflection * Lt / Lc);
31       NotDone = Standard_False;
32     }
33   }
34 }
35
36
37 //=======================================================================
38 //function : GCPnts_TangentialDeflection
39 //purpose  : 
40 //=======================================================================
41
42 GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
43  const TheCurve&        C,
44  const Standard_Real    AngularDeflection,
45  const Standard_Real    CurvatureDeflection,
46  const Standard_Integer MinimumOfPoints,
47  const Standard_Real    UTol)
48
49
50   Initialize (C,AngularDeflection,CurvatureDeflection,MinimumOfPoints,UTol); 
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    FirstParameter,
62  const Standard_Real    LastParameter,
63  const Standard_Real    AngularDeflection,
64  const Standard_Real    CurvatureDeflection,
65  const Standard_Integer MinimumOfPoints,
66  const Standard_Real    UTol)
67
68
69   Initialize (C, 
70               FirstParameter, 
71               LastParameter,
72               AngularDeflection, 
73               CurvatureDeflection, 
74               MinimumOfPoints,
75               UTol);
76 }
77
78
79
80 //=======================================================================
81 //function : Initialize
82 //purpose  : 
83 //=======================================================================
84
85 void GCPnts_TangentialDeflection::Initialize (
86  const TheCurve&        C, 
87  const Standard_Real    AngularDeflection, 
88  const Standard_Real    CurvatureDeflection,
89  const Standard_Integer MinimumOfPoints,
90  const Standard_Real    UTol)
91
92 {
93   Initialize (C, 
94               C.FirstParameter (), 
95               C.LastParameter (),
96               AngularDeflection, 
97               CurvatureDeflection,
98               MinimumOfPoints,
99               UTol);
100 }
101
102
103 //=======================================================================
104 //function : Initialize
105 //purpose  : 
106 //=======================================================================
107
108 void GCPnts_TangentialDeflection::Initialize (
109  const TheCurve&        C, 
110  const Standard_Real    FirstParameter, 
111  const Standard_Real    LastParameter,
112  const Standard_Real    AngularDeflection, 
113  const Standard_Real    CurvatureDeflection,
114  const Standard_Integer MinimumOfPoints,
115  const Standard_Real    UTol)
116
117 {
118   
119   Standard_ConstructionError_Raise_if (CurvatureDeflection <= Precision::Confusion () || AngularDeflection   <= Precision::Angular (), "GCPnts_TangentialDeflection::Initialize - Zero Deflection")
120
121  parameters.Clear ();
122  points    .Clear ();
123  if (FirstParameter < LastParameter) {
124    firstu = FirstParameter;
125    lastu  = LastParameter;
126  }
127  else {
128    lastu  = FirstParameter;
129    firstu = LastParameter;
130  }
131  uTol                = UTol;
132  angularDeflection   = AngularDeflection;
133  curvatureDeflection = CurvatureDeflection;
134  minNbPnts           = Max (MinimumOfPoints, 2);
135
136  switch (C.GetType()) {
137
138  case GeomAbs_Line:   
139    PerformLinear (C);
140    break;
141
142  case GeomAbs_Circle: 
143    PerformCircular (C);
144    break;
145
146  case GeomAbs_BSplineCurve:
147    {
148      Handle_TheBSplineCurve BS = C.BSpline() ;
149      if (BS->NbPoles() == 2 ) PerformLinear (C);
150      else                     PerformCurve  (C);
151      break;
152    }
153  case GeomAbs_BezierCurve:
154    {
155      Handle_TheBezierCurve  BZ = C.Bezier(); 
156      if (BZ->NbPoles() == 2) PerformLinear (C);
157      else                    PerformCurve  (C);
158      break;
159    }
160  default: PerformCurve (C);
161    
162  }
163 }
164
165
166 //=======================================================================
167 //function : PerformLinear
168 //purpose  : 
169 //=======================================================================
170
171 void GCPnts_TangentialDeflection::PerformLinear (const TheCurve& C) {
172
173   gp_Pnt P;
174   D0 (C, firstu, P);
175   parameters.Append (firstu);
176   points    .Append (P);
177   if (minNbPnts > 2) {
178     Standard_Real Du = (lastu - firstu) / minNbPnts;
179     Standard_Real U = firstu + Du;
180     for (Standard_Integer i = 2; i <= minNbPnts; i++) {
181       D0 (C, U, P);
182       parameters.Append (U);
183       points    .Append (P);
184       U += Du;
185     }
186   }
187   D0 (C, lastu, P);
188   parameters.Append (lastu);
189   points    .Append (P);
190 }
191
192
193 //=======================================================================
194 //function : PerformCircular
195 //purpose  : 
196 //=======================================================================
197
198 void GCPnts_TangentialDeflection::PerformCircular (const TheCurve& C) 
199 {
200   // akm 8/01/02 : check the radius before divide by it
201   Standard_Real dfR = C.Circle().Radius();
202   Standard_Real Du = 0.;
203   if (Abs(dfR) > Precision::Confusion())
204     Du = Max(1.0e0 - (curvatureDeflection/dfR),0.0e0) ;
205   Du  = acos (Du); Du+=Du;
206   Du               = Min (Du, angularDeflection);
207   Standard_Integer NbPoints = (Standard_Integer )((lastu - firstu) / Du);
208   NbPoints         = Max (NbPoints, minNbPnts-1);
209   Du               = (lastu - firstu) / NbPoints;
210
211   gp_Pnt P;
212   Standard_Real U = firstu;
213   for (Standard_Integer i = 1; i <= NbPoints; i++) {
214     D0 (C, U, P);
215     parameters.Append (U);
216     points    .Append (P);
217     U += Du;
218   }
219   D0 (C, lastu, P);
220   parameters.Append (lastu);
221   points    .Append (P);
222 }
223
224
225
226 //=======================================================================
227 //function : PerformCurve
228 //purpose  : On respecte ll'angle et la fleche, on peut imposer un nombre
229 //           minimum de points sur un element lineaire
230 //=======================================================================
231
232 void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
233
234 {
235   Standard_Integer i;       
236   Standard_Boolean Correction, TooLarge, TooSmall, MorePoints;
237   gp_XYZ V1, V2;
238 //  gp_Pnt MiddlePoint, CurrentPoint, LastPoint, PreviousPoint;   
239   gp_Pnt MiddlePoint, CurrentPoint, LastPoint;   
240   Standard_Real Coef, ACoef, FCoef, Du, Dusave, MiddleU, L1, L2;
241
242   Standard_Real U1       = firstu;   
243   Standard_Real U2       = firstu;   
244   Standard_Real AngleMax = angularDeflection * 0.5;  //car on prend le point milieu
245   Standard_Real LTol     = Precision::Confusion ();  //protection longueur nulle
246   Standard_Real ATol     = Precision::Angular ();    //protection angle nul
247
248   D0 (C, lastu, LastPoint);
249
250   //Initialization du calcul
251   //---------------------------------------------------------------------
252
253   Standard_Boolean NotDone = Standard_True;
254   Dusave = (lastu - firstu)*Us3;
255   Du     = Dusave;
256   EvaluateDu (C, U1, CurrentPoint, Du, NotDone);
257   parameters.Append (U1);
258   points    .Append (CurrentPoint);
259
260   if (NotDone) {
261     //C'est soit une droite, soit une singularite :
262     V1 = LastPoint.XYZ ();
263     V1.Subtract (CurrentPoint.XYZ());
264     L1 = V1.Modulus ();
265     if (L1 > LTol) {
266       //Si c'est une droite on verifie en calculant minNbPoints :
267       Standard_Boolean IsLine   = Standard_True;
268       Standard_Integer NbPoints = 3;
269       if (minNbPnts > 3) NbPoints = minNbPnts;
270       Du = (lastu-firstu)/NbPoints;
271       MiddleU = firstu + Du;
272       for (i = 2; i < NbPoints; i++) {
273         D0 (C, MiddleU, MiddlePoint);
274         V2 = MiddlePoint.XYZ();
275         V2.Subtract (CurrentPoint.XYZ());
276         L2 = V2.Modulus ();
277         if (L2 > LTol) {
278           if (((V2.CrossMagnitude (V1))/(L1*L2)) >= ATol) {
279             //C'etait une singularite
280             IsLine = Standard_False;
281             break;
282           }
283           if (minNbPnts > 2) {
284             parameters.Append (MiddleU);
285             points    .Append (MiddlePoint);
286           }
287         }
288         MiddleU += Du;
289       }
290       if (IsLine) {
291         //C'etait une droite (plusieurs poles alignes), Calcul termine :
292         parameters.Append (lastu);
293         points    .Append (LastPoint);
294         return;
295       }
296       else {
297         //c'etait une singularite on continue :
298         Standard_Integer pointsLength=points.Length ();
299         for (i = 2; i <= pointsLength; i++) {
300           points    .Remove (i);
301           parameters.Remove (i);
302           pointsLength--;
303         }
304         Du = (lastu-firstu)*Us3;
305       }
306     }
307     else  {
308       
309       Du = (lastu-firstu)/2.1;
310       MiddleU = firstu + Du;
311       D0 (C, MiddleU, MiddlePoint);
312       V1 = MiddlePoint.XYZ ();
313       V1.Subtract (CurrentPoint.XYZ());
314       L1 = V1.Modulus ();
315       if (L1 < LTol) {
316         // L1 < LTol C'est une courbe de longueur nulle, calcul termine :
317         // on renvoi un segment de 2 points   (protection)
318         parameters.Append (lastu);
319         points    .Append (LastPoint);
320         return;
321       }
322     }
323   }
324
325   if (Du > Dusave) Du = Dusave;
326   else             Dusave = Du;
327
328   if (Du < uTol) {
329     Du = lastu - firstu;
330     if (Du < uTol) {
331       parameters.Append (lastu);
332       points    .Append (LastPoint);
333       return;
334     }
335   }
336
337   //Traitement normal pour une courbe
338   //-------------------------------------------------------------------
339   MorePoints = Standard_True;
340
341   while (MorePoints) {
342
343     U2 += Du;                             
344
345     if (U2 >= lastu) {                       //Bout de courbe
346       U2 = lastu;
347       CurrentPoint = LastPoint;
348       Du = U2-U1;
349       Dusave = Du;
350     }
351     else D0 (C, U2, CurrentPoint);           //Point suivant
352
353     TooLarge   = Standard_False;
354     TooSmall   = Standard_False;
355     Correction = Standard_True;
356      
357     while (Correction) {                     //Ajustement Du
358       MiddleU = (U1+U2)*0.5;                 //Verif / au point milieu
359       D0 (C, MiddleU, MiddlePoint);
360
361       V1 = CurrentPoint.XYZ ();              //Critere de fleche
362       V1.Subtract (points.Last().XYZ());
363       V2 = MiddlePoint.XYZ ();
364       V2.Subtract (points.Last().XYZ());
365       L1 = V1.Modulus ();
366       if (L1 > LTol) FCoef = V1.CrossMagnitude(V2)/(L1*curvatureDeflection);
367       else           FCoef = 0.0;
368
369       V1 = CurrentPoint.XYZ ();              //Critere d'angle
370       V1.Subtract (MiddlePoint.XYZ ());
371       L1 = V1.Modulus ();
372       L2 = V2.Modulus ();
373       Standard_Real angg = V1.CrossMagnitude(V2)/(L1*L2);
374       if (L1 > LTol && L2 > LTol) ACoef = angg/AngleMax;
375       else                        ACoef = 0.0;
376
377       if (ACoef >= FCoef) Coef = ACoef;      //On retient le plus penalisant
378       else                Coef = FCoef;
379
380
381       if (Coef <= 1.0) {                     
382         if (Abs (lastu-U2) < uTol) {
383           parameters.Append (lastu);
384           points    .Append (LastPoint);
385           MorePoints = Standard_False;
386           Correction = Standard_False;
387         }
388         else {
389           if (Coef >= 0.75 || TooLarge || TooSmall) { 
390             parameters.Append (U2);
391             points    .Append (CurrentPoint);
392             Correction = Standard_False;
393           }
394           else {
395             TooSmall = Standard_True;
396             Standard_Real UUU2 = U2;
397             Du += Min((U2-U1)*(1.-Coef), Du*Us3);
398
399             U2 = U1 + Du;
400             if (U2 >= lastu) U2 = UUU2;
401             if (U2 >= lastu) {
402               parameters.Append (lastu);
403               points    .Append (LastPoint);
404               MorePoints = Standard_False;
405               Correction = Standard_False;
406             }
407             else D0 (C, U2, CurrentPoint);
408           }
409         }
410       }
411       else {
412
413         if (Coef >= 1.5) {
414           U2 = MiddleU;
415           CurrentPoint = MiddlePoint;
416         }
417         else {
418           Du*=0.9;
419           U2 = U1 + Du;
420           D0 (C, U2, CurrentPoint);
421           TooLarge = Standard_True;
422         }
423
424       }
425     }
426     
427     Du  = U2-U1;
428
429     if (MorePoints) {
430       if (U1 > firstu) {
431         if (FCoef > ACoef) {
432           //La fleche est critere de decoupage
433           EvaluateDu (C, U2, CurrentPoint, Du, NotDone);
434           if (NotDone) {
435             Du += (Du-Dusave)*(Du/Dusave);
436             if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
437             if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
438           }
439         }
440         else {
441           //L'angle est le critere de decoupage
442           Du += (Du-Dusave)*(Du/Dusave);
443           if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
444           if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
445         }
446       }
447       
448       if (Du < uTol) {
449         Du = lastu - U2;
450         if (Du < uTol) {
451           parameters.Append (lastu);
452           points    .Append (LastPoint);
453           MorePoints = Standard_False;
454         }
455         else if (Du*Us3 > uTol) Du*=Us3;
456       }
457       U1 = U2;
458       Dusave = Du;
459     }
460   }
461                               //Recalage avant dernier point :
462   i = points.Length()-1;
463 //  Real d = points (i).Distance (points (i+1));
464 // if (Abs(parameters (i) - parameters (i+1))<= 0.000001 || d < Precision::Confusion()) {
465 //    cout<<"deux points confondus"<<endl;
466 //    parameters.Remove (i+1);
467 //    points.Remove (i+1);
468 //    i--;
469 //  }
470   
471   if (i >= 2) {
472     MiddleU = parameters (i-1);
473     MiddleU = (lastu + MiddleU)*0.5;
474     D0 (C, MiddleU, MiddlePoint);
475     parameters.SetValue (i, MiddleU);
476     points    .SetValue (i, MiddlePoint);
477   }
478
479   //-- ------------------------------------------------------------
480   //-- On rajoute des points aux milieux des segments si le nombre
481   //-- mini de points n'est pas atteint
482   //--
483   Standard_Integer Nbp =  points.Length();
484   Standard_Integer MinNb= (9*minNbPnts)/10;
485   if(MinNb<3) MinNb=3;
486
487   //-- if(Nbp <  MinNb) { cout<<"\n*"; } else {  cout<<"\n."; } 
488   while(Nbp < MinNb) { 
489     //-- cout<<" \nGCPnts TangentialDeflection : Ajout de Points ("<<Nbp<<" "<<minNbPnts<<" )"<<endl;
490     for(i=2; i<=Nbp; i++) { 
491       MiddleU = (parameters.Value(i-1)+parameters.Value(i))*0.5;
492       D0 (C, MiddleU, MiddlePoint); 
493       parameters.InsertBefore(i,MiddleU);
494       points.InsertBefore(i,MiddlePoint);
495       Nbp++;
496       i++;
497     }
498   }
499 }
500
501
502
503
504