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