0026106: BRepMesh - revision of data model
[occt.git] / src / GCPnts / GCPnts_TangentialDeflection.pxx
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 #include <NCollection_List.hxx>
25 #include <math_PSO.hxx>
26 #include <math_BrentMinimum.hxx>
27
28 #define Us3 0.3333333333333333333333333333
29
30 void GCPnts_TangentialDeflection::EvaluateDu (
31   const TheCurve& C,
32   const Standard_Real      U,
33   gp_Pnt&   P,
34   Standard_Real&     Du,
35   Standard_Boolean&  NotDone) const {
36
37   gp_Vec T, N;
38   D2 (C, U, P, T, N);
39   Standard_Real Lt   = T.Magnitude ();
40   Standard_Real LTol = Precision::Confusion ();
41   if (Lt > LTol && N.Magnitude () > LTol) {
42     Standard_Real Lc = N.CrossMagnitude (T);
43     Standard_Real Ln = Lc/Lt;
44     if (Ln > LTol) {
45       Du = sqrt (8.0 * Max(curvatureDeflection, myMinLen) / Ln);
46       NotDone = Standard_False;
47     }
48   }
49 }
50
51
52 //=======================================================================
53 //function : GCPnts_TangentialDeflection
54 //purpose  : 
55 //=======================================================================
56
57 GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
58  const TheCurve&        C,
59  const Standard_Real    AngularDeflection,
60  const Standard_Real    CurvatureDeflection,
61  const Standard_Integer MinimumOfPoints,
62  const Standard_Real    UTol,
63  const Standard_Real    theMinLen)
64
65
66   Initialize (C,AngularDeflection,CurvatureDeflection,MinimumOfPoints,UTol,theMinLen); 
67 }
68
69
70 //=======================================================================
71 //function : GCPnts_TangentialDeflection
72 //purpose  : 
73 //=======================================================================
74
75 GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
76  const TheCurve&        C,
77  const Standard_Real    FirstParameter,
78  const Standard_Real    LastParameter,
79  const Standard_Real    AngularDeflection,
80  const Standard_Real    CurvatureDeflection,
81  const Standard_Integer MinimumOfPoints,
82  const Standard_Real    UTol,
83  const Standard_Real    theMinLen)
84
85
86   Initialize (C, 
87         FirstParameter, 
88         LastParameter,
89         AngularDeflection, 
90         CurvatureDeflection, 
91         MinimumOfPoints,
92         UTol, theMinLen);
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  const Standard_Real    theMinLen)
109
110 {
111   Initialize (C, 
112               C.FirstParameter (), 
113               C.LastParameter (),
114               AngularDeflection, 
115               CurvatureDeflection,
116               MinimumOfPoints,
117               UTol, theMinLen);
118 }
119
120
121 //=======================================================================
122 //function : Initialize
123 //purpose  : 
124 //=======================================================================
125
126 void GCPnts_TangentialDeflection::Initialize (
127  const TheCurve&        C, 
128  const Standard_Real    FirstParameter, 
129  const Standard_Real    LastParameter,
130  const Standard_Real    AngularDeflection, 
131  const Standard_Real    CurvatureDeflection,
132  const Standard_Integer MinimumOfPoints,
133  const Standard_Real    UTol,
134  const Standard_Real    theMinLen)
135
136 {
137   
138   Standard_ConstructionError_Raise_if (CurvatureDeflection <= Precision::Confusion () || AngularDeflection   <= Precision::Angular (), "GCPnts_TangentialDeflection::Initialize - Zero Deflection")
139
140  parameters.Clear ();
141  points    .Clear ();
142  if (FirstParameter < LastParameter) {
143    firstu = FirstParameter;
144    lastu  = LastParameter;
145  }
146  else {
147    lastu  = FirstParameter;
148    firstu = LastParameter;
149  }
150  uTol                = UTol;
151  angularDeflection   = AngularDeflection;
152  curvatureDeflection = CurvatureDeflection;
153  minNbPnts           = Max (MinimumOfPoints, 2);
154  myMinLen            = Max(theMinLen, Precision::Confusion());
155
156  switch (C.GetType()) {
157
158  case GeomAbs_Line:   
159    PerformLinear (C);
160    break;
161
162  case GeomAbs_Circle: 
163    PerformCircular (C);
164    break;
165
166  case GeomAbs_BSplineCurve:
167    {
168      Handle_TheBSplineCurve BS = C.BSpline() ;
169      if (BS->NbPoles() == 2 ) PerformLinear (C);
170      else                     PerformCurve  (C);
171      break;
172    }
173  case GeomAbs_BezierCurve:
174    {
175      Handle_TheBezierCurve  BZ = C.Bezier(); 
176      if (BZ->NbPoles() == 2) PerformLinear (C);
177      else                    PerformCurve  (C);
178      break;
179    }
180  default: PerformCurve (C);
181    
182  }
183 }
184
185
186 //=======================================================================
187 //function : PerformLinear
188 //purpose  : 
189 //=======================================================================
190
191 void GCPnts_TangentialDeflection::PerformLinear (const TheCurve& C) {
192
193   gp_Pnt P;
194   D0 (C, firstu, P);
195   parameters.Append (firstu);
196   points    .Append (P);
197   if (minNbPnts > 2) {
198     Standard_Real Du = (lastu - firstu) / minNbPnts;
199     Standard_Real U = firstu + Du;
200     for (Standard_Integer i = 2; i <= minNbPnts; i++) {
201       D0 (C, U, P);
202       parameters.Append (U);
203       points    .Append (P);
204       U += Du;
205     }
206   }
207   D0 (C, lastu, P);
208   parameters.Append (lastu);
209   points    .Append (P);
210 }
211
212 //=======================================================================
213 //function : PerformCircular
214 //purpose  : 
215 //=======================================================================
216
217 void GCPnts_TangentialDeflection::PerformCircular (const TheCurve& C) 
218 {
219   // akm 8/01/02 : check the radius before divide by it
220   Standard_Real dfR = C.Circle().Radius();
221   Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep(
222     dfR, curvatureDeflection, angularDeflection, myMinLen);
223     
224   const Standard_Real aDiff = lastu - firstu;
225   // Round up number of points to satisfy curvatureDeflection more precisely
226   Standard_Integer NbPoints = (Standard_Integer)Min(Ceiling(aDiff / Du), 1.0e+6);
227   NbPoints = Max(NbPoints, minNbPnts - 1);
228   Du       = aDiff / NbPoints;
229
230   gp_Pnt P;
231   Standard_Real U = firstu;
232   for (Standard_Integer i = 1; i <= NbPoints; i++)
233   {
234     D0 (C, U, P);
235     parameters.Append (U);
236     points    .Append (P);
237     U += Du;
238   }
239   D0 (C, lastu, P);
240   parameters.Append (lastu);
241   points    .Append (P);
242 }
243
244
245 //=======================================================================
246 //function : PerformCurve
247 //purpose  : On respecte ll'angle et la fleche, on peut imposer un nombre
248 //           minimum de points sur un element lineaire
249 //=======================================================================
250 void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
251
252 {
253   Standard_Integer i, j;       
254   gp_XYZ V1, V2;
255   gp_Pnt MiddlePoint, CurrentPoint, LastPoint;   
256   Standard_Real Du, Dusave, MiddleU, L1, L2;
257
258   Standard_Real U1       = firstu;   
259   Standard_Real LTol     = Precision::Confusion ();  //protection longueur nulle
260   Standard_Real ATol     = 1.e-2 * angularDeflection;
261   if(ATol > 1.e-2)
262     ATol = 1.e-2;
263   else if(ATol < 1.e-7)
264     ATol = 1.e-7;
265
266   D0 (C, lastu, LastPoint);
267
268   //Initialization du calcul
269
270   Standard_Boolean NotDone = Standard_True;
271   Dusave = (lastu - firstu)*Us3;
272   Du     = Dusave;
273   EvaluateDu (C, U1, CurrentPoint, Du, NotDone);
274   parameters.Append (U1);
275   points    .Append (CurrentPoint);
276
277   // Used to detect "isLine" current bspline and in Du computation in general handling.
278   Standard_Integer NbInterv = C.NbIntervals(GeomAbs_CN);
279   TColStd_Array1OfReal Intervs(1, NbInterv+1);
280   C.Intervals(Intervs, GeomAbs_CN);
281
282   if (NotDone || Du > 5. * Dusave) {
283     //C'est soit une droite, soit une singularite :
284     V1 = (LastPoint.XYZ() - CurrentPoint.XYZ());
285     L1 = V1.Modulus ();
286     if (L1 > LTol)
287     {
288       //Si c'est une droite on verifie en calculant minNbPoints :
289       Standard_Boolean IsLine   = Standard_True;
290       Standard_Integer NbPoints = (minNbPnts > 3) ? minNbPnts : 3;
291       switch (C.GetType()) {
292       case GeomAbs_BSplineCurve:
293         {
294           Handle_TheBSplineCurve BS = C.BSpline() ;
295           NbPoints = Max(BS->Degree() + 1, NbPoints);
296           break;
297         }
298       case GeomAbs_BezierCurve:
299         {
300           Handle_TheBezierCurve  BZ = C.Bezier();
301           NbPoints = Max(BZ->Degree() + 1, NbPoints);
302           break;
303         }
304       default:
305       ;}
306       ////
307       Standard_Real param = 0.;
308       for (i = 1; i <= NbInterv && IsLine; ++i)
309       {
310         // Avoid usage intervals out of [firstu, lastu].
311         if ((Intervs(i+1) < firstu) ||
312             (Intervs(i)   > lastu))
313         {
314           continue;
315         }
316         // Fix border points in applicable intervals, to avoid be out of target interval.
317         if ((Intervs(i)   < firstu) &&
318             (Intervs(i+1) > firstu))
319         {
320           Intervs(i) = firstu;
321         }
322         if ((Intervs(i)   < lastu) &&
323             (Intervs(i+1) > lastu))
324         {
325           Intervs(i + 1) = lastu;
326         }
327
328         Standard_Real delta = (Intervs(i+1) - Intervs(i))/NbPoints;
329         for (j = 1; j <= NbPoints && IsLine; ++j)
330         {
331           param = Intervs(i) + j*delta;
332           D0 (C, param, MiddlePoint);
333           V2 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
334           L2 = V2.Modulus ();
335           if (L2 > LTol)
336           {
337             const Standard_Real aAngle = V2.CrossMagnitude(V1)/(L1*L2);
338             IsLine = (aAngle < ATol);
339           }
340         }
341       }
342
343       if (IsLine)
344       {
345         parameters.Clear();
346         points    .Clear();
347
348         PerformLinear(C);
349         return;
350       }
351       else
352       {
353         //c'etait une singularite on continue :
354         //Du = Dusave;
355         EvaluateDu (C, param, MiddlePoint, Du, NotDone);
356       }
357     }
358     else
359     {
360       Du = (lastu-firstu)/2.1;
361       MiddleU = firstu + Du;
362       D0 (C, MiddleU, MiddlePoint);
363       V1 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
364       L1 = V1.Modulus ();
365       if (L1 < LTol)
366       {
367         // L1 < LTol C'est une courbe de longueur nulle, calcul termine :
368         // on renvoi un segment de 2 points   (protection)
369         parameters.Append (lastu);
370         points    .Append (LastPoint);
371         return;
372       }
373     }
374   }
375
376   if (Du > Dusave) Du = Dusave;
377   else             Dusave = Du;
378
379   if (Du < uTol) {
380     Du = lastu - firstu;
381     if (Du < uTol) {
382       parameters.Append (lastu);
383       points    .Append (LastPoint);
384       return;
385     }
386   }
387
388   //Traitement normal pour une courbe
389   Standard_Boolean MorePoints = Standard_True;
390   Standard_Real U2            = firstu;   
391   Standard_Real AngleMax      = angularDeflection * 0.5;  //car on prend le point milieu
392   Standard_Integer aIdx[2] = {Intervs.Lower(), Intervs.Lower()}; // Indexes of intervals of U1 and U2, used to handle non-uniform case.
393   Standard_Boolean isNeedToCheck = Standard_False;
394   gp_Pnt aPrevPoint = points.Last();
395
396   while (MorePoints) {
397     aIdx[0] = getIntervalIdx(U1, Intervs, aIdx[0]);
398     U2 += Du;
399
400     if (U2 >= lastu) {                       //Bout de courbe
401       U2 = lastu;
402       CurrentPoint = LastPoint;
403       Du = U2-U1;
404       Dusave = Du;
405     }
406     else D0 (C, U2, CurrentPoint);           //Point suivant
407
408     Standard_Real Coef = 0.0, ACoef = 0., FCoef = 0.;
409     Standard_Boolean Correction, TooLarge, TooSmall;
410     TooLarge   = Standard_False;
411     Correction = Standard_True;
412     TooSmall = Standard_False;
413
414     while (Correction) {                     //Ajustement Du
415       if (isNeedToCheck)
416       {
417         aIdx[1] = getIntervalIdx(U2, Intervs, aIdx[0]);
418         if (aIdx[1] > aIdx[0]) // Jump to another polynom.
419         {
420           if (Du > (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3) // Set Du to the smallest value and check deflection on it.
421           {
422             Du = (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3;
423             U2 = U1 + Du;
424             if (U2 > lastu)
425               U2 = lastu;
426             D0 (C, U2, CurrentPoint);
427           }
428         }
429       }
430       MiddleU = (U1+U2)*0.5;                 //Verif / au point milieu
431       D0 (C, MiddleU, MiddlePoint);
432
433       V1 = (CurrentPoint.XYZ() - aPrevPoint.XYZ()); //Critere de fleche
434       V2 = (MiddlePoint.XYZ()  - aPrevPoint.XYZ());
435       L1 = V1.Modulus ();
436
437       FCoef = (L1 > myMinLen) ? 
438         V1.CrossMagnitude(V2)/(L1*curvatureDeflection) : 0.0;
439
440       V1 = (CurrentPoint.XYZ() - MiddlePoint.XYZ()); //Critere d'angle
441       L1 = V1.Modulus ();
442       L2 = V2.Modulus ();
443       if (L1 > myMinLen && L2 > myMinLen)
444       {
445         Standard_Real angg = V1.CrossMagnitude(V2) / (L1 * L2);
446         ACoef = angg / AngleMax;
447       }
448       else
449         ACoef = 0.0;
450
451       //On retient le plus penalisant
452       Coef = Max(ACoef, FCoef);
453
454       if (isNeedToCheck && Coef < 0.55)
455       {
456         isNeedToCheck = Standard_False;
457         Du = Dusave;
458         U2 = U1 + Du;
459         if (U2 > lastu)
460           U2 = lastu;
461         D0 (C, U2, CurrentPoint);
462         continue;
463       }
464
465       if (Coef <= 1.0) {
466         if (Abs (lastu-U2) < uTol) {
467           parameters.Append (lastu);
468           points    .Append (LastPoint);
469           MorePoints = Standard_False;
470           Correction = Standard_False;
471         }
472         else {
473           if (Coef >= 0.55 || TooLarge) { 
474             parameters.Append (U2);
475             points    .Append (CurrentPoint);
476             aPrevPoint = CurrentPoint;
477             Correction = Standard_False;
478             isNeedToCheck = Standard_True;
479           }
480           else if (TooSmall) {
481             Correction = Standard_False;
482             aPrevPoint = CurrentPoint;
483           }
484           else {
485             TooSmall = Standard_True;
486             //Standard_Real UUU2 = U2;
487             Du += Min((U2-U1)*(1.-Coef), Du*Us3);
488
489             U2 = U1 + Du;
490             if (U2 > lastu)
491               U2 = lastu;
492             D0 (C, U2, CurrentPoint);
493           }
494         }
495       }
496       else {
497
498         if (Coef >= 1.5) {
499           if (!aPrevPoint.IsEqual(points.Last(), Precision::Confusion()))
500           {
501             parameters.Append (U1);
502             points    .Append (aPrevPoint);
503           }
504           U2 = MiddleU;
505           Du  = U2-U1;
506           CurrentPoint = MiddlePoint;
507         }
508         else {
509           Du*=0.9;
510           U2 = U1 + Du;
511           D0 (C, U2, CurrentPoint);
512           TooLarge = Standard_True;
513         }
514
515       }
516     }
517     
518     Du  = U2-U1;
519
520     if (MorePoints) {
521       if (U1 > firstu) {
522         if (FCoef > ACoef) {
523           //La fleche est critere de decoupage
524           EvaluateDu (C, U2, CurrentPoint, Du, NotDone);
525           if (NotDone) {
526             Du += (Du-Dusave)*(Du/Dusave);
527             if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
528             if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
529           }
530         }
531         else {
532           //L'angle est le critere de decoupage
533           Du += (Du-Dusave)*(Du/Dusave);
534           if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
535           if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
536         }
537       }
538       
539       if (Du < uTol) {
540         Du = lastu - U2;
541         if (Du < uTol) {
542           parameters.Append (lastu);
543           points    .Append (LastPoint);
544           MorePoints = Standard_False;
545         }
546         else if (Du*Us3 > uTol) Du*=Us3;
547       }
548       U1 = U2;
549       Dusave = Du;
550     }
551   }
552   //Recalage avant dernier point :
553   i = points.Length()-1;
554 //  Real d = points (i).Distance (points (i+1));
555 // if (Abs(parameters (i) - parameters (i+1))<= 0.000001 || d < Precision::Confusion()) {
556 //    cout<<"deux points confondus"<<endl;
557 //    parameters.Remove (i+1);
558 //    points.Remove (i+1);
559 //    i--;
560 //  }
561   if (i >= 2) {
562     MiddleU = parameters (i-1);
563     MiddleU = (lastu + MiddleU)*0.5;
564     D0 (C, MiddleU, MiddlePoint);
565     parameters.SetValue (i, MiddleU);
566     points    .SetValue (i, MiddlePoint);
567   }
568
569   //-- On rajoute des points aux milieux des segments si le nombre
570   //-- mini de points n'est pas atteint
571   //--
572   Standard_Integer Nbp =  points.Length();
573   Standard_Integer MinNb= (9*minNbPnts)/10;
574   //if(MinNb<4) MinNb=4;
575
576   //-- if(Nbp <  MinNb) { cout<<"\n*"; } else {  cout<<"\n."; } 
577   while(Nbp < MinNb) { 
578     //-- cout<<" \nGCPnts TangentialDeflection : Ajout de Points ("<<Nbp<<" "<<minNbPnts<<" )"<<endl;
579     for (i = 2; i <= Nbp; i += 2) {
580       MiddleU = (parameters.Value(i-1)+parameters.Value(i))*0.5;
581       D0 (C, MiddleU, MiddlePoint); 
582       parameters.InsertBefore(i,MiddleU);
583       points.InsertBefore(i,MiddlePoint);
584       Nbp++;
585     }
586   }
587   //Additional check for intervals
588   Standard_Real MinLen2 = myMinLen * myMinLen;
589   Standard_Integer MaxNbp = 10 * Nbp;
590   for(i = 1; i < Nbp; ++i)
591   {
592     U1 = parameters(i);
593     U2 = parameters(i + 1);
594     // Check maximal deflection on interval;
595     Standard_Real dmax = 0.;
596     Standard_Real umax = 0.;
597     Standard_Real amax = 0.;
598     EstimDefl(C, U1, U2, dmax, umax);
599     const gp_Pnt& P1 = points(i);
600     const gp_Pnt& P2 = points(i+1);
601     D0(C, umax, MiddlePoint);
602     amax = EstimAngl(P1, MiddlePoint, P2);
603     if(dmax > curvatureDeflection || amax > AngleMax)
604     {
605       if(umax - U1 > uTol && U2 - umax > uTol)
606       {
607         if (P1.SquareDistance(MiddlePoint) > MinLen2 && P2.SquareDistance(MiddlePoint) > MinLen2)
608         {
609           parameters.InsertAfter(i, umax);
610           points.InsertAfter(i, MiddlePoint);
611           ++Nbp;
612           --i; //To compensate ++i in loop header: i must point to first part of splitted interval
613           if(Nbp > MaxNbp)
614           {
615             break;
616           }
617         }
618       }
619     }
620   }
621   // 
622 }
623
624 //=======================================================================
625 //function : EstimDefl
626 //purpose  : Estimation of maximal deflection for interval [U1, U2]
627 //           
628 //=======================================================================
629 void GCPnts_TangentialDeflection::EstimDefl (const TheCurve& C,
630                             const Standard_Real U1, const Standard_Real U2, 
631                             Standard_Real& MaxDefl, Standard_Real& UMax)
632 {
633   Standard_Real Du = (lastu - firstu);
634   //
635   TheMaxCurvLinDist aFunc(C, U1, U2);
636   //
637   const Standard_Integer aNbIter = 100;
638   Standard_Real reltol = Max(1.e-3, 2.*uTol/((Abs(U1) + Abs(U2))));
639   //
640   math_BrentMinimum anOptLoc(reltol, aNbIter, uTol);
641   anOptLoc.Perform(aFunc, U1, (U1+U2)/2., U2);
642   if(anOptLoc.IsDone())
643   {
644     MaxDefl = Sqrt(-anOptLoc.Minimum());
645     UMax = anOptLoc.Location();
646     return;
647   }
648   //
649   math_Vector aLowBorder(1,1);
650   math_Vector aUppBorder(1,1);
651   math_Vector aSteps(1,1);
652   //
653   aSteps(1) = Max(0.1 * Du, 100. * uTol);
654   Standard_Integer aNbParticles = Max(8, RealToInt(32 * (U2 - U1) / Du));
655   //
656   aLowBorder(1) = U1;
657   aUppBorder(1) = U2;
658   //
659   //
660   Standard_Real aValue;
661   math_Vector aT(1,1);
662   TheMaxCurvLinDistMV aFuncMV(aFunc);
663
664   math_PSO aFinder(&aFuncMV, aLowBorder, aUppBorder, aSteps, aNbParticles); 
665   aFinder.Perform(aSteps, aValue, aT);
666   //
667   anOptLoc.Perform(aFunc, Max(aT(1) - aSteps(1), U1) , aT(1), Min(aT(1) + aSteps(1),U2));
668   if(anOptLoc.IsDone())
669   {
670     MaxDefl = Sqrt(-anOptLoc.Minimum());
671     UMax = anOptLoc.Location();
672     return;
673   }
674   MaxDefl = Sqrt(-aValue);
675   UMax = aT(1);
676   //
677 }