0a657c07c21c4d133086c646752f74f7f7953279
[occt.git] / src / GCPnts / GCPnts_AbscissaPoint.gxx
1 // Created on: 1995-05-05
2 // Created by: Modelistation
3 // Copyright (c) 1995-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 // Dimension independant used to implement GCPnts_AbscissaPoint
18
19 // compute the type 
20 // and the length ratio if GCPnts_LengthParametrized
21 #include <GCPnts_AbscissaType.hxx>
22 #include <gp_Vec.hxx>
23 #include <gp_Vec2d.hxx>
24 #include <gp_Circ.hxx>
25 #include <gp_Circ2d.hxx>
26 #include <Precision.hxx>
27 #include <TColStd_Array1OfReal.hxx>
28 #include <BSplCLib.hxx>
29
30 static GCPnts_AbscissaType computeType( TheCurve& C,
31                                        Standard_Real& Ratio)
32 {
33   GCPnts_AbscissaType LocalType ;
34
35   if (C.NbIntervals(GeomAbs_CN) > 1)
36     return GCPnts_AbsComposite;
37
38   switch (C.GetType()) {
39     
40   case GeomAbs_Line:
41     Ratio = 1.0e0 ;
42     return GCPnts_LengthParametrized;
43
44   case GeomAbs_Circle:
45     Ratio = C.Circle().Radius();
46     return GCPnts_LengthParametrized;
47
48   case GeomAbs_BezierCurve:
49     {
50       Handle_TheBezierCurve Bz = C.Bezier();
51       if ((Bz->NbPoles() == 2) && !(Bz->IsRational())) {
52         Ratio = Bz->DN(0,1).Magnitude();
53         LocalType =  GCPnts_LengthParametrized;
54       }
55       else
56         LocalType = GCPnts_Parametrized;
57       return LocalType ;
58     }
59   case GeomAbs_BSplineCurve:
60     {
61       Handle_TheBSplineCurve Bs = C.BSpline();
62       if ((Bs->NbPoles() == 2) && !(Bs->IsRational())) {
63         Ratio = Bs->DN(Bs->FirstParameter(),1).Magnitude();
64         LocalType = GCPnts_LengthParametrized;
65       }
66       else
67         LocalType = GCPnts_Parametrized;
68       return LocalType ; 
69     }
70     default:
71       return GCPnts_Parametrized;
72     
73   }
74 }
75
76 // compute a point at distance Abscis from parameter U0
77 // using Ui as initial guess
78
79 static void Compute(CPnts_AbscissaPoint& theComputer,
80                     TheCurve& C,
81                     Standard_Real& Abscis,
82                     Standard_Real& U0,
83                     Standard_Real& Ui,
84                     const Standard_Real EPSILON) 
85 {
86   // test for easy solution
87   if (Abs(Abscis) <= Precision::Confusion()) {
88     theComputer.SetParameter(U0);
89     return;
90   }
91
92   Standard_Real Ratio;
93   GCPnts_AbscissaType Type = computeType(C,Ratio);
94
95   switch (Type) {
96   case GCPnts_LengthParametrized : 
97     theComputer.SetParameter(U0 + Abscis / Ratio);
98     return;
99
100   case GCPnts_Parametrized : 
101     theComputer.Init(C);
102     theComputer.Perform(Abscis, U0, Ui, EPSILON);
103     return;
104
105   case GCPnts_AbsComposite : 
106     {
107       Standard_Integer NbIntervals = C.NbIntervals(GeomAbs_CN);
108       TColStd_Array1OfReal TI(1,NbIntervals+1);
109       C.Intervals(TI,GeomAbs_CN);
110       Standard_Real L = 0.0, sign = 1.;
111       Standard_Integer Index = 1;
112       BSplCLib::Hunt(TI,U0,Index);
113       Standard_Integer Direction = 1;
114       if (Abscis < 0) {
115         Direction = 0;
116         Abscis = -Abscis;
117         sign = -1.;
118       }
119
120       while ((Index >= 1) && (Index <= NbIntervals)) {
121
122         L = CPnts_AbscissaPoint::Length(C, U0, TI(Index+Direction));
123         if (Abs(L - Abscis) <= Precision::Confusion()) {
124           theComputer.SetParameter(TI(Index+Direction));
125           return;
126         }
127         if(L > Abscis) {
128           if ((Ui < TI(Index)) || (Ui > TI(Index+1))) {
129             Ui = (Abscis / L) * (TI(Index+1) - U0);
130             if (Direction)
131               Ui = U0 + Ui;
132             else
133               Ui = U0 - Ui;
134           }
135           theComputer.Init(C,TI(Index),TI(Index+1));
136           theComputer.Perform(sign*Abscis, U0, Ui, EPSILON);
137           return;
138         }
139         else {
140           U0 = TI(Index+Direction);
141           Abscis -= L;
142         }
143         if (Direction)
144           Index++;
145         else
146           Index--;
147       }
148
149       // Push a little bit outside the limits (hairy !!!)
150       Ui = U0 + 0.1;
151       theComputer.Init(C,U0,U0+0.2);
152       theComputer.Perform(sign*Abscis, U0, Ui, EPSILON);
153       return;
154     }
155     break;
156   }
157
158 }
159
160 // introduced by rbv for curvilinear parametrization
161 // performs more apropriate tolerance managment
162
163 static void AdvCompute(CPnts_AbscissaPoint& theComputer,
164                     TheCurve& C,
165                     Standard_Real& Abscis,
166                     Standard_Real& U0,
167                     Standard_Real& Ui,
168                     const Standard_Real EPSILON) 
169 {
170   // test for easy solution
171   if (Abs(Abscis) <= EPSILON) {
172     theComputer.SetParameter(U0);
173     return;
174   }
175   
176   Standard_Real Ratio;
177   GCPnts_AbscissaType Type = computeType(C,Ratio);
178
179   switch (Type) {
180   case GCPnts_LengthParametrized : 
181     theComputer.SetParameter(U0 + Abscis / Ratio);
182     return;
183
184   case GCPnts_Parametrized : 
185 //    theComputer.Init(C);
186     theComputer.Init(C, EPSILON); //rbv's modification
187 //
188     theComputer.AdvPerform(Abscis, U0, Ui, EPSILON);
189     return;
190
191   case GCPnts_AbsComposite : 
192     {
193       Standard_Integer NbIntervals = C.NbIntervals(GeomAbs_CN);
194       TColStd_Array1OfReal TI(1,NbIntervals+1);
195       C.Intervals(TI,GeomAbs_CN);
196       Standard_Real L = 0.0, sign = 1.;
197       Standard_Integer Index = 1;
198       BSplCLib::Hunt(TI,U0,Index);
199         
200       Standard_Integer Direction = 1;
201       if (Abscis < 0) {
202         Direction = 0;
203         Abscis = -Abscis;
204         sign = -1.;
205       }
206
207       if(Index == 0 && Direction > 0) {
208         L = CPnts_AbscissaPoint::Length(C, U0, TI(Index+Direction), EPSILON);
209         if (Abs(L - Abscis) <= /*Precision::Confusion()*/EPSILON) {
210           theComputer.SetParameter(TI(Index+Direction));
211           return;
212         }
213         if(L > Abscis) {
214           if ( Ui > TI(Index+1) ) {
215             Ui = (Abscis / L) * (TI(Index+1) - U0);
216             Ui = U0 + Ui;
217           }
218           theComputer.Init(C,U0,TI(Index+1), EPSILON);
219           theComputer.AdvPerform(sign*Abscis, U0, Ui, EPSILON);
220           return;
221         }
222         else {
223           U0 = TI(Index+Direction);
224           Abscis -= L;
225         }
226         Index++;
227       }
228           
229
230       while ((Index >= 1) && (Index <= NbIntervals)) {
231
232         L = CPnts_AbscissaPoint::Length(C, U0, TI(Index+Direction), EPSILON);
233         if (Abs(L - Abscis) <= /*Precision::Confusion()*/EPSILON) {
234           theComputer.SetParameter(TI(Index+Direction));
235           return;
236         }
237         if(L > Abscis) {
238           if ((Ui < TI(Index)) || (Ui > TI(Index+1))) {
239             Ui = (Abscis / L) * (TI(Index+1) - U0);
240             if (Direction)
241               Ui = U0 + Ui;
242             else
243               Ui = U0 - Ui;
244           }
245           theComputer.Init(C,TI(Index),TI(Index+1), EPSILON);
246           theComputer.AdvPerform(sign*Abscis, U0, Ui, EPSILON);
247           return;
248         }
249         else {
250           U0 = TI(Index+Direction);
251           Abscis -= L;
252         }
253         if (Direction) {
254           Index++;
255
256         }
257         else {
258           Index--;
259
260         }
261       }
262
263       // Push a little bit outside the limits (hairy !!!)
264
265       Standard_Boolean nonperiodic = !C.IsPeriodic();
266       Ui = U0 + sign*0.1;
267       Standard_Real U1 = U0 + sign*.2;
268       if(nonperiodic) {
269         if(sign > 0) {
270           Ui = Min(Ui,C.LastParameter());
271           U1 = Min(U1, C.LastParameter());
272         }
273         else {
274           Ui = Max(Ui,C.FirstParameter());
275           U1 = Max(U1, C.FirstParameter());
276         }
277       }
278           
279       theComputer.Init(C, U0, U1, EPSILON);
280       theComputer.AdvPerform(sign*Abscis, U0, Ui, EPSILON);
281       return;
282     }
283     break;
284   }
285
286 }
287
288 //=======================================================================
289 //function : Length
290 //purpose  : 
291 //=======================================================================
292
293 Standard_Real GCPnts_AbscissaPoint::Length(TheCurve& C)
294 {
295   return GCPnts_AbscissaPoint::Length(C,C.FirstParameter(), 
296                                       C.LastParameter());
297 }
298
299 //=======================================================================
300 //function : Length
301 //purpose  : 
302 //=======================================================================
303
304 Standard_Real GCPnts_AbscissaPoint::Length(TheCurve& C,
305                                            const Standard_Real Tol)
306 {
307   return GCPnts_AbscissaPoint::Length(C,C.FirstParameter(), 
308                                       C.LastParameter(),Tol);
309 }
310
311
312 //=======================================================================
313 //function : Length
314 //purpose  : 
315 //=======================================================================
316
317 Standard_Real GCPnts_AbscissaPoint::Length(TheCurve& C,
318                                            const Standard_Real U1,
319                                            const Standard_Real U2)
320 {
321   Standard_Real Ratio;
322   GCPnts_AbscissaType Type = computeType(C,Ratio);
323   switch (Type) {
324
325   case GCPnts_LengthParametrized: 
326     return Abs(U2-U1) * Ratio;
327                 
328   case GCPnts_Parametrized: 
329     return CPnts_AbscissaPoint::Length(C, U1, U2);
330
331   case GCPnts_AbsComposite: 
332     {
333       Standard_Integer NbIntervals = C.NbIntervals(GeomAbs_CN);
334       TColStd_Array1OfReal TI(1,NbIntervals+1);
335       C.Intervals(TI,GeomAbs_CN);
336       Standard_Real UU1 = Min(U1, U2);
337       Standard_Real UU2 = Max(U1, U2);
338       Standard_Real L = 0.0;
339       for(Standard_Integer Index = 1; Index <= NbIntervals; Index++) {
340         if (TI(Index)   > UU2) break;
341         if (TI(Index+1) < UU1) continue;
342         L += CPnts_AbscissaPoint::Length(C, 
343                                          Max(TI(Index),UU1),
344                                          Min(TI(Index+1),UU2));
345       }
346       return L;
347     }
348   }
349   return RealLast();
350 }
351
352 //=======================================================================
353 //function : Length
354 //purpose  : 
355 //=======================================================================
356
357 Standard_Real GCPnts_AbscissaPoint::Length(TheCurve& C,
358                                            const Standard_Real U1,
359                                            const Standard_Real U2,
360                                            const Standard_Real Tol)
361 {
362   Standard_Real Ratio;
363   GCPnts_AbscissaType Type = computeType(C,Ratio);
364   switch (Type) {
365
366   case GCPnts_LengthParametrized: 
367     return Abs(U2-U1) * Ratio;
368                 
369   case GCPnts_Parametrized: 
370     return CPnts_AbscissaPoint::Length(C, U1, U2, Tol);
371
372   case GCPnts_AbsComposite: 
373     {
374       Standard_Integer NbIntervals = C.NbIntervals(GeomAbs_CN);
375       TColStd_Array1OfReal TI(1,NbIntervals+1);
376       C.Intervals(TI,GeomAbs_CN);
377       Standard_Real UU1 = Min(U1, U2);
378       Standard_Real UU2 = Max(U1, U2);
379       Standard_Real L = 0.0;
380       for(Standard_Integer Index = 1; Index <= NbIntervals; Index++) {
381         if (TI(Index)   > UU2) break;
382         if (TI(Index+1) < UU1) continue;
383         L += CPnts_AbscissaPoint::Length(C, 
384                                          Max(TI(Index),UU1),
385                                          Min(TI(Index+1),UU2),
386                                          Tol);
387       }
388       return L;
389     }
390   }
391   return RealLast();
392 }
393
394
395 //=======================================================================
396 //function : GCPnts_AbscissaPoint
397 //purpose  : 
398 //=======================================================================
399
400 GCPnts_AbscissaPoint::GCPnts_AbscissaPoint
401                                   (TheCurve& C,
402                                    const Standard_Real   Abscissa,
403                                    const Standard_Real   U0)
404 {
405   Standard_Real L = GCPnts_AbscissaPoint::Length(C);
406   if (L < Precision::Confusion()) {
407     Standard_ConstructionError::Raise();
408   } 
409   Standard_Real Abscis = Abscissa;
410   Standard_Real UU0 = U0;
411   Standard_Real UUi = U0 + 
412     (Abscis / L) * (C.LastParameter() - C.FirstParameter());
413   Compute(myComputer, C, Abscis, UU0, UUi,
414           C.Resolution(Precision::Confusion()));  
415 }
416
417 //=======================================================================
418 //function : GCPnts_AbscissaPoint
419 //purpose  : rbv for curvilinear parametrization
420 //=======================================================================
421
422 GCPnts_AbscissaPoint::GCPnts_AbscissaPoint
423                                   (const Standard_Real Tol,
424                                    TheCurve& C,
425                                    const Standard_Real   Abscissa,
426                                    const Standard_Real   U0)
427 {
428   Standard_Real L = GCPnts_AbscissaPoint::Length(C, Tol);
429 /*  if (L < Precision::Confusion()) {
430     cout<<"FirstParameter = "<<C.FirstParameter()<<endl;
431     cout<<"LastParameter = "<<C.LastParameter()<<endl;
432     Standard_ConstructionError::Raise("GCPnts_AbscissaPoint::GCPnts_AbscissaPoint");
433   } 
434 */
435   Standard_Real Abscis = Abscissa;
436   Standard_Real UU0 = U0;
437   Standard_Real UUi;
438   if (L >= Precision::Confusion()) 
439     UUi= U0 + 
440       (Abscis / L) * (C.LastParameter() - C.FirstParameter());
441   else UUi = U0;
442
443   AdvCompute(myComputer, C, Abscis, UU0, UUi, Tol);  
444 }
445
446 //=======================================================================
447 //function : GCPnts_AbscissaPoint
448 //purpose  : 
449 //=======================================================================
450
451 GCPnts_AbscissaPoint::GCPnts_AbscissaPoint
452                                   (TheCurve& C,
453                                    const Standard_Real   Abscissa,
454                                    const Standard_Real   U0,
455                                    const Standard_Real   Ui)
456 {
457   Standard_Real Abscis = Abscissa;
458   Standard_Real UU0 = U0;
459   Standard_Real UUi = Ui;
460   Compute(myComputer, C, Abscis, UU0, UUi, 
461           C.Resolution(Precision::Confusion()));
462 }
463
464 //=======================================================================
465 //function : GCPnts_AbscissaPoint
466 //purpose  : rbv for curvilinear parametrization
467 //=======================================================================
468
469 GCPnts_AbscissaPoint::GCPnts_AbscissaPoint
470                                   (TheCurve& C,
471                                    const Standard_Real   Abscissa,
472                                    const Standard_Real   U0,
473                                    const Standard_Real   Ui,
474                                    const Standard_Real   Tol)
475 {
476   Standard_Real Abscis = Abscissa;
477   Standard_Real UU0 = U0;
478   Standard_Real UUi = Ui;
479   AdvCompute(myComputer, C, Abscis, UU0, UUi, Tol);
480 }