017e43f3dc3a2bbe2eff2a04e778e8b4e78800c0
[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 under
9 // the terms of the GNU Lesser General Public License 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 * Max(curvatureDeflection, myMinLen) / 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  const Standard_Real    theMinLen)
61
62
63   Initialize (C,AngularDeflection,CurvatureDeflection,MinimumOfPoints,UTol,theMinLen); 
64 }
65
66
67 //=======================================================================
68 //function : GCPnts_TangentialDeflection
69 //purpose  : 
70 //=======================================================================
71
72 GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
73  const TheCurve&        C,
74  const Standard_Real    FirstParameter,
75  const Standard_Real    LastParameter,
76  const Standard_Real    AngularDeflection,
77  const Standard_Real    CurvatureDeflection,
78  const Standard_Integer MinimumOfPoints,
79  const Standard_Real    UTol,
80  const Standard_Real    theMinLen)
81
82
83   Initialize (C, 
84         FirstParameter, 
85         LastParameter,
86         AngularDeflection, 
87         CurvatureDeflection, 
88         MinimumOfPoints,
89         UTol, theMinLen);
90 }
91
92
93
94 //=======================================================================
95 //function : Initialize
96 //purpose  : 
97 //=======================================================================
98
99 void GCPnts_TangentialDeflection::Initialize (
100  const TheCurve&        C, 
101  const Standard_Real    AngularDeflection, 
102  const Standard_Real    CurvatureDeflection,
103  const Standard_Integer MinimumOfPoints,
104  const Standard_Real    UTol,
105  const Standard_Real    theMinLen)
106
107 {
108   Initialize (C, 
109               C.FirstParameter (), 
110               C.LastParameter (),
111               AngularDeflection, 
112               CurvatureDeflection,
113               MinimumOfPoints,
114               UTol, theMinLen);
115 }
116
117
118 //=======================================================================
119 //function : Initialize
120 //purpose  : 
121 //=======================================================================
122
123 void GCPnts_TangentialDeflection::Initialize (
124  const TheCurve&        C, 
125  const Standard_Real    FirstParameter, 
126  const Standard_Real    LastParameter,
127  const Standard_Real    AngularDeflection, 
128  const Standard_Real    CurvatureDeflection,
129  const Standard_Integer MinimumOfPoints,
130  const Standard_Real    UTol,
131  const Standard_Real    theMinLen)
132
133 {
134   
135   Standard_ConstructionError_Raise_if (CurvatureDeflection <= Precision::Confusion () || AngularDeflection   <= Precision::Angular (), "GCPnts_TangentialDeflection::Initialize - Zero Deflection")
136
137  parameters.Clear ();
138  points    .Clear ();
139  if (FirstParameter < LastParameter) {
140    firstu = FirstParameter;
141    lastu  = LastParameter;
142  }
143  else {
144    lastu  = FirstParameter;
145    firstu = LastParameter;
146  }
147  uTol                = UTol;
148  angularDeflection   = AngularDeflection;
149  curvatureDeflection = CurvatureDeflection;
150  minNbPnts           = Max (MinimumOfPoints, 2);
151  myMinLen            = Max(theMinLen, Precision::Confusion());
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 //function : PerformCircular
211 //purpose  : 
212 //=======================================================================
213
214 void GCPnts_TangentialDeflection::PerformCircular (const TheCurve& C) 
215 {
216   // akm 8/01/02 : check the radius before divide by it
217   Standard_Real dfR = C.Circle().Radius();
218   Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep(
219     dfR, curvatureDeflection, angularDeflection, myMinLen);
220     
221   const Standard_Real aDiff = lastu - firstu;
222   Standard_Integer NbPoints = (Standard_Integer)(aDiff / Du);
223   NbPoints = Max(NbPoints, minNbPnts - 1);
224   Du       = aDiff / NbPoints;
225
226   gp_Pnt P;
227   Standard_Real U = firstu;
228   for (Standard_Integer i = 1; i <= NbPoints; i++)
229   {
230     D0 (C, U, P);
231     parameters.Append (U);
232     points    .Append (P);
233     U += Du;
234   }
235   D0 (C, lastu, P);
236   parameters.Append (lastu);
237   points    .Append (P);
238 }
239
240
241 //=======================================================================
242 //function : PerformCurve
243 //purpose  : On respecte ll'angle et la fleche, on peut imposer un nombre
244 //           minimum de points sur un element lineaire
245 //=======================================================================
246 void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
247
248 {
249   Standard_Integer i, j;       
250   gp_XYZ V1, V2;
251   gp_Pnt MiddlePoint, CurrentPoint, LastPoint;   
252   Standard_Real Du, Dusave, MiddleU, L1, L2;
253
254   Standard_Real U1       = firstu;   
255   Standard_Real LTol     = Precision::Confusion ();  //protection longueur nulle
256   Standard_Real ATol     = Precision::Angular ();    //protection angle nul
257
258   D0 (C, lastu, LastPoint);
259
260   //Initialization du calcul
261
262   Standard_Boolean NotDone = Standard_True;
263   Dusave = (lastu - firstu)*Us3;
264   Du     = Dusave;
265   EvaluateDu (C, U1, CurrentPoint, Du, NotDone);
266   parameters.Append (U1);
267   points    .Append (CurrentPoint);
268
269   // Used to detect "isLine" current bspline and in Du computation in general handling.
270   Standard_Integer NbInterv = C.NbIntervals(GeomAbs_CN);
271   TColStd_Array1OfReal Intervs(1, NbInterv+1);
272   C.Intervals(Intervs, GeomAbs_CN);
273
274   if (NotDone) {
275     //C'est soit une droite, soit une singularite :
276     V1 = (LastPoint.XYZ() - CurrentPoint.XYZ());
277     L1 = V1.Modulus ();
278     if (L1 > LTol)
279     {
280       //Si c'est une droite on verifie en calculant minNbPoints :
281       Standard_Boolean IsLine   = Standard_True;
282       Standard_Integer NbPoints = (minNbPnts > 3) ? minNbPnts : 3;
283       ////
284       Standard_Real param = 0.;
285       for (i = 1; i <= NbInterv && IsLine; ++i)
286       {
287         // Avoid usage intervals out of [firstu, lastu].
288         if ((Intervs(i+1) < firstu) ||
289             (Intervs(i)   > lastu))
290         {
291           continue;
292         }
293         // Fix border points in applicable intervals, to avoid be out of target interval.
294         if ((Intervs(i)   < firstu) &&
295             (Intervs(i+1) > firstu))
296         {
297           Intervs(i) = firstu;
298         }
299         if ((Intervs(i)   < lastu) &&
300             (Intervs(i+1) > lastu))
301         {
302           Intervs(i + 1) = lastu;
303         }
304
305         Standard_Real delta = (Intervs(i+1) - Intervs(i))/NbPoints;
306         for (j = 1; j <= NbPoints && IsLine; ++j)
307         {
308           param = Intervs(i) + j*delta;
309           D0 (C, param, MiddlePoint);
310           V2 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
311           L2 = V2.Modulus ();
312           if (L2 > LTol)
313           {
314             const Standard_Real aAngle = V2.CrossMagnitude(V1)/(L1*L2);
315             IsLine = (aAngle < ATol);
316           }
317         }
318       }
319
320       if (IsLine)
321       {
322         parameters.Clear();
323         points    .Clear();
324
325         PerformLinear(C);
326         return;
327       }
328       else
329       {
330         //c'etait une singularite on continue :
331         //Du = Dusave;
332         EvaluateDu (C, param, MiddlePoint, Du, NotDone);
333       }
334     }
335     else
336     {
337       Du = (lastu-firstu)/2.1;
338       MiddleU = firstu + Du;
339       D0 (C, MiddleU, MiddlePoint);
340       V1 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
341       L1 = V1.Modulus ();
342       if (L1 < LTol)
343       {
344         // L1 < LTol C'est une courbe de longueur nulle, calcul termine :
345         // on renvoi un segment de 2 points   (protection)
346         parameters.Append (lastu);
347         points    .Append (LastPoint);
348         return;
349       }
350     }
351   }
352
353   if (Du > Dusave) Du = Dusave;
354   else             Dusave = Du;
355
356   if (Du < uTol) {
357     Du = lastu - firstu;
358     if (Du < uTol) {
359       parameters.Append (lastu);
360       points    .Append (LastPoint);
361       return;
362     }
363   }
364
365   //Traitement normal pour une courbe
366   Standard_Boolean MorePoints = Standard_True;
367   Standard_Real U2            = firstu;   
368   Standard_Real AngleMax      = angularDeflection * 0.5;  //car on prend le point milieu
369   Standard_Integer aIdx[2] = {Intervs.Lower(), Intervs.Lower()}; // Indexes of intervals of U1 and U2, used to handle non-uniform case.
370   Standard_Boolean isNeedToCheck = Standard_False;
371   gp_Pnt aPrevPoint = points.Last();
372
373   while (MorePoints) {
374     aIdx[0] = getIntervalIdx(U1, Intervs, aIdx[0]);
375     U2 += Du;
376
377     if (U2 >= lastu) {                       //Bout de courbe
378       U2 = lastu;
379       CurrentPoint = LastPoint;
380       Du = U2-U1;
381       Dusave = Du;
382     }
383     else D0 (C, U2, CurrentPoint);           //Point suivant
384
385     Standard_Real Coef = 0.0, ACoef = 0., FCoef = 0.;
386     Standard_Boolean Correction, TooLarge, TooSmall;
387     TooLarge   = Standard_False;
388     Correction = Standard_True;
389     TooSmall = Standard_False;
390
391     while (Correction) {                     //Ajustement Du
392       if (isNeedToCheck)
393       {
394         aIdx[1] = getIntervalIdx(U2, Intervs, aIdx[0]);
395         if (aIdx[1] > aIdx[0]) // Jump to another polynom.
396         {
397           if (Du > (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3) // Set Du to the smallest value and check deflection on it.
398           {
399             Du = (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3;
400             U2 = U1 + Du;
401             if (U2 > lastu)
402               U2 = lastu;
403             D0 (C, U2, CurrentPoint);
404           }
405         }
406       }
407       MiddleU = (U1+U2)*0.5;                 //Verif / au point milieu
408       D0 (C, MiddleU, MiddlePoint);
409
410       V1 = (CurrentPoint.XYZ() - aPrevPoint.XYZ()); //Critere de fleche
411       V2 = (MiddlePoint.XYZ()  - aPrevPoint.XYZ());
412       L1 = V1.Modulus ();
413
414       FCoef = (L1 > myMinLen) ? 
415         V1.CrossMagnitude(V2)/(L1*curvatureDeflection) : 0.0;
416
417       V1 = (CurrentPoint.XYZ() - MiddlePoint.XYZ()); //Critere d'angle
418       L1 = V1.Modulus ();
419       L2 = V2.Modulus ();
420       if (L1 > myMinLen && L2 > myMinLen)
421       {
422         Standard_Real angg = V1.CrossMagnitude(V2) / (L1 * L2);
423         ACoef = angg / AngleMax;
424       }
425       else
426         ACoef = 0.0;
427
428       //On retient le plus penalisant
429       Coef = Max(ACoef, FCoef);
430
431       if (isNeedToCheck && Coef < 0.55)
432       {
433         isNeedToCheck = Standard_False;
434         Du = Dusave;
435         U2 = U1 + Du;
436         if (U2 > lastu)
437           U2 = lastu;
438         D0 (C, U2, CurrentPoint);
439         continue;
440       }
441
442       if (Coef <= 1.0) {
443         if (Abs (lastu-U2) < uTol) {
444           parameters.Append (lastu);
445           points    .Append (LastPoint);
446           MorePoints = Standard_False;
447           Correction = Standard_False;
448         }
449         else {
450           if (Coef >= 0.55 || TooLarge) { 
451             parameters.Append (U2);
452             points    .Append (CurrentPoint);
453             aPrevPoint = CurrentPoint;
454             Correction = Standard_False;
455             isNeedToCheck = Standard_True;
456           }
457           else if (TooSmall) {
458             Correction = Standard_False;
459             aPrevPoint = CurrentPoint;
460           }
461           else {
462             TooSmall = Standard_True;
463             //Standard_Real UUU2 = U2;
464             Du += Min((U2-U1)*(1.-Coef), Du*Us3);
465
466             U2 = U1 + Du;
467             if (U2 > lastu)
468               U2 = lastu;
469             D0 (C, U2, CurrentPoint);
470           }
471         }
472       }
473       else {
474
475         if (Coef >= 1.5) {
476           U2 = MiddleU;
477           Du  = U2-U1;
478           CurrentPoint = MiddlePoint;
479         }
480         else {
481           Du*=0.9;
482           U2 = U1 + Du;
483           D0 (C, U2, CurrentPoint);
484           TooLarge = Standard_True;
485         }
486
487       }
488     }
489     
490     Du  = U2-U1;
491
492     if (MorePoints) {
493       if (U1 > firstu) {
494         if (FCoef > ACoef) {
495           //La fleche est critere de decoupage
496           EvaluateDu (C, U2, CurrentPoint, Du, NotDone);
497           if (NotDone) {
498             Du += (Du-Dusave)*(Du/Dusave);
499             if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
500             if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
501           }
502         }
503         else {
504           //L'angle est le critere de decoupage
505           Du += (Du-Dusave)*(Du/Dusave);
506           if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
507           if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
508         }
509       }
510       
511       if (Du < uTol) {
512         Du = lastu - U2;
513         if (Du < uTol) {
514           parameters.Append (lastu);
515           points    .Append (LastPoint);
516           MorePoints = Standard_False;
517         }
518         else if (Du*Us3 > uTol) Du*=Us3;
519       }
520       U1 = U2;
521       Dusave = Du;
522     }
523   }
524   //Recalage avant dernier point :
525   i = points.Length()-1;
526 //  Real d = points (i).Distance (points (i+1));
527 // if (Abs(parameters (i) - parameters (i+1))<= 0.000001 || d < Precision::Confusion()) {
528 //    cout<<"deux points confondus"<<endl;
529 //    parameters.Remove (i+1);
530 //    points.Remove (i+1);
531 //    i--;
532 //  }
533   if (i >= 2) {
534     MiddleU = parameters (i-1);
535     MiddleU = (lastu + MiddleU)*0.5;
536     D0 (C, MiddleU, MiddlePoint);
537     parameters.SetValue (i, MiddleU);
538     points    .SetValue (i, MiddlePoint);
539   }
540
541   //-- On rajoute des points aux milieux des segments si le nombre
542   //-- mini de points n'est pas atteint
543   //--
544   Standard_Integer Nbp =  points.Length();
545   Standard_Integer MinNb= (9*minNbPnts)/10;
546   //if(MinNb<4) MinNb=4;
547
548   //-- if(Nbp <  MinNb) { cout<<"\n*"; } else {  cout<<"\n."; } 
549   while(Nbp < MinNb) { 
550     //-- cout<<" \nGCPnts TangentialDeflection : Ajout de Points ("<<Nbp<<" "<<minNbPnts<<" )"<<endl;
551     for(i=2; i<=Nbp; i++) { 
552       MiddleU = (parameters.Value(i-1)+parameters.Value(i))*0.5;
553       D0 (C, MiddleU, MiddlePoint); 
554       parameters.InsertBefore(i,MiddleU);
555       points.InsertBefore(i,MiddlePoint);
556       Nbp++;
557       i++;
558     }
559   }
560 }