0024474: GCPnts_AbscissaPoint calculates invalid point
[occt.git] / src / GCPnts / GCPnts_AbscissaPoint.gxx
CommitLineData
b311480e 1// Created on: 1995-05-05
2// Created by: Modelistation
3// Copyright (c) 1995-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
7fd59977 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
973c2be1 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.
7fd59977 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
b311480e 16
17// Dimension independant used to implement GCPnts_AbscissaPoint
7fd59977 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
30static 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
79static 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
163static 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{
7fd59977 170 Standard_Real Ratio;
171 GCPnts_AbscissaType Type = computeType(C,Ratio);
172
173 switch (Type) {
174 case GCPnts_LengthParametrized :
175 theComputer.SetParameter(U0 + Abscis / Ratio);
176 return;
177
178 case GCPnts_Parametrized :
179// theComputer.Init(C);
180 theComputer.Init(C, EPSILON); //rbv's modification
181//
182 theComputer.AdvPerform(Abscis, U0, Ui, EPSILON);
183 return;
184
185 case GCPnts_AbsComposite :
186 {
187 Standard_Integer NbIntervals = C.NbIntervals(GeomAbs_CN);
188 TColStd_Array1OfReal TI(1,NbIntervals+1);
189 C.Intervals(TI,GeomAbs_CN);
190 Standard_Real L = 0.0, sign = 1.;
191 Standard_Integer Index = 1;
192 BSplCLib::Hunt(TI,U0,Index);
193
194 Standard_Integer Direction = 1;
195 if (Abscis < 0) {
196 Direction = 0;
197 Abscis = -Abscis;
198 sign = -1.;
199 }
200
201 if(Index == 0 && Direction > 0) {
202 L = CPnts_AbscissaPoint::Length(C, U0, TI(Index+Direction), EPSILON);
203 if (Abs(L - Abscis) <= /*Precision::Confusion()*/EPSILON) {
204 theComputer.SetParameter(TI(Index+Direction));
205 return;
206 }
207 if(L > Abscis) {
208 if ( Ui > TI(Index+1) ) {
209 Ui = (Abscis / L) * (TI(Index+1) - U0);
210 Ui = U0 + Ui;
211 }
212 theComputer.Init(C,U0,TI(Index+1), EPSILON);
213 theComputer.AdvPerform(sign*Abscis, U0, Ui, EPSILON);
214 return;
215 }
216 else {
217 U0 = TI(Index+Direction);
218 Abscis -= L;
219 }
220 Index++;
221 }
222
223
224 while ((Index >= 1) && (Index <= NbIntervals)) {
225
226 L = CPnts_AbscissaPoint::Length(C, U0, TI(Index+Direction), EPSILON);
bb0e6b9b 227 if (Abs(L - Abscis) <= Precision::PConfusion()) {
7fd59977 228 theComputer.SetParameter(TI(Index+Direction));
229 return;
230 }
231 if(L > Abscis) {
232 if ((Ui < TI(Index)) || (Ui > TI(Index+1))) {
233 Ui = (Abscis / L) * (TI(Index+1) - U0);
234 if (Direction)
235 Ui = U0 + Ui;
236 else
237 Ui = U0 - Ui;
238 }
239 theComputer.Init(C,TI(Index),TI(Index+1), EPSILON);
240 theComputer.AdvPerform(sign*Abscis, U0, Ui, EPSILON);
241 return;
242 }
243 else {
244 U0 = TI(Index+Direction);
245 Abscis -= L;
246 }
247 if (Direction) {
248 Index++;
249
250 }
251 else {
252 Index--;
253
254 }
255 }
256
257 // Push a little bit outside the limits (hairy !!!)
258
259 Standard_Boolean nonperiodic = !C.IsPeriodic();
260 Ui = U0 + sign*0.1;
261 Standard_Real U1 = U0 + sign*.2;
262 if(nonperiodic) {
263 if(sign > 0) {
264 Ui = Min(Ui,C.LastParameter());
265 U1 = Min(U1, C.LastParameter());
266 }
267 else {
268 Ui = Max(Ui,C.FirstParameter());
269 U1 = Max(U1, C.FirstParameter());
270 }
271 }
272
273 theComputer.Init(C, U0, U1, EPSILON);
274 theComputer.AdvPerform(sign*Abscis, U0, Ui, EPSILON);
275 return;
276 }
277 break;
278 }
279
280}
281
282//=======================================================================
283//function : Length
284//purpose :
285//=======================================================================
286
287Standard_Real GCPnts_AbscissaPoint::Length(TheCurve& C)
288{
289 return GCPnts_AbscissaPoint::Length(C,C.FirstParameter(),
290 C.LastParameter());
291}
292
293//=======================================================================
294//function : Length
295//purpose :
296//=======================================================================
297
298Standard_Real GCPnts_AbscissaPoint::Length(TheCurve& C,
299 const Standard_Real Tol)
300{
301 return GCPnts_AbscissaPoint::Length(C,C.FirstParameter(),
302 C.LastParameter(),Tol);
303}
304
305
306//=======================================================================
307//function : Length
308//purpose :
309//=======================================================================
310
311Standard_Real GCPnts_AbscissaPoint::Length(TheCurve& C,
312 const Standard_Real U1,
313 const Standard_Real U2)
314{
315 Standard_Real Ratio;
316 GCPnts_AbscissaType Type = computeType(C,Ratio);
317 switch (Type) {
318
319 case GCPnts_LengthParametrized:
320 return Abs(U2-U1) * Ratio;
321
322 case GCPnts_Parametrized:
323 return CPnts_AbscissaPoint::Length(C, U1, U2);
324
325 case GCPnts_AbsComposite:
326 {
327 Standard_Integer NbIntervals = C.NbIntervals(GeomAbs_CN);
328 TColStd_Array1OfReal TI(1,NbIntervals+1);
329 C.Intervals(TI,GeomAbs_CN);
330 Standard_Real UU1 = Min(U1, U2);
331 Standard_Real UU2 = Max(U1, U2);
332 Standard_Real L = 0.0;
333 for(Standard_Integer Index = 1; Index <= NbIntervals; Index++) {
334 if (TI(Index) > UU2) break;
335 if (TI(Index+1) < UU1) continue;
336 L += CPnts_AbscissaPoint::Length(C,
337 Max(TI(Index),UU1),
338 Min(TI(Index+1),UU2));
339 }
340 return L;
341 }
342 }
343 return RealLast();
344}
345
346//=======================================================================
347//function : Length
348//purpose :
349//=======================================================================
350
351Standard_Real GCPnts_AbscissaPoint::Length(TheCurve& C,
352 const Standard_Real U1,
353 const Standard_Real U2,
354 const Standard_Real Tol)
355{
356 Standard_Real Ratio;
357 GCPnts_AbscissaType Type = computeType(C,Ratio);
358 switch (Type) {
359
360 case GCPnts_LengthParametrized:
361 return Abs(U2-U1) * Ratio;
362
363 case GCPnts_Parametrized:
364 return CPnts_AbscissaPoint::Length(C, U1, U2, Tol);
365
366 case GCPnts_AbsComposite:
367 {
368 Standard_Integer NbIntervals = C.NbIntervals(GeomAbs_CN);
369 TColStd_Array1OfReal TI(1,NbIntervals+1);
370 C.Intervals(TI,GeomAbs_CN);
371 Standard_Real UU1 = Min(U1, U2);
372 Standard_Real UU2 = Max(U1, U2);
373 Standard_Real L = 0.0;
374 for(Standard_Integer Index = 1; Index <= NbIntervals; Index++) {
375 if (TI(Index) > UU2) break;
376 if (TI(Index+1) < UU1) continue;
377 L += CPnts_AbscissaPoint::Length(C,
378 Max(TI(Index),UU1),
379 Min(TI(Index+1),UU2),
380 Tol);
381 }
382 return L;
383 }
384 }
385 return RealLast();
386}
387
388
389//=======================================================================
390//function : GCPnts_AbscissaPoint
391//purpose :
392//=======================================================================
393
394GCPnts_AbscissaPoint::GCPnts_AbscissaPoint
395 (TheCurve& C,
396 const Standard_Real Abscissa,
397 const Standard_Real U0)
398{
399 Standard_Real L = GCPnts_AbscissaPoint::Length(C);
400 if (L < Precision::Confusion()) {
401 Standard_ConstructionError::Raise();
402 }
403 Standard_Real Abscis = Abscissa;
404 Standard_Real UU0 = U0;
405 Standard_Real UUi = U0 +
406 (Abscis / L) * (C.LastParameter() - C.FirstParameter());
407 Compute(myComputer, C, Abscis, UU0, UUi,
408 C.Resolution(Precision::Confusion()));
409}
410
411//=======================================================================
412//function : GCPnts_AbscissaPoint
413//purpose : rbv for curvilinear parametrization
414//=======================================================================
415
416GCPnts_AbscissaPoint::GCPnts_AbscissaPoint
417 (const Standard_Real Tol,
418 TheCurve& C,
419 const Standard_Real Abscissa,
420 const Standard_Real U0)
421{
422 Standard_Real L = GCPnts_AbscissaPoint::Length(C, Tol);
423/* if (L < Precision::Confusion()) {
424 cout<<"FirstParameter = "<<C.FirstParameter()<<endl;
425 cout<<"LastParameter = "<<C.LastParameter()<<endl;
426 Standard_ConstructionError::Raise("GCPnts_AbscissaPoint::GCPnts_AbscissaPoint");
427 }
428*/
429 Standard_Real Abscis = Abscissa;
430 Standard_Real UU0 = U0;
431 Standard_Real UUi;
432 if (L >= Precision::Confusion())
433 UUi= U0 +
434 (Abscis / L) * (C.LastParameter() - C.FirstParameter());
435 else UUi = U0;
436
437 AdvCompute(myComputer, C, Abscis, UU0, UUi, Tol);
438}
439
440//=======================================================================
441//function : GCPnts_AbscissaPoint
442//purpose :
443//=======================================================================
444
445GCPnts_AbscissaPoint::GCPnts_AbscissaPoint
446 (TheCurve& C,
447 const Standard_Real Abscissa,
448 const Standard_Real U0,
449 const Standard_Real Ui)
450{
451 Standard_Real Abscis = Abscissa;
452 Standard_Real UU0 = U0;
453 Standard_Real UUi = Ui;
454 Compute(myComputer, C, Abscis, UU0, UUi,
455 C.Resolution(Precision::Confusion()));
456}
457
458//=======================================================================
459//function : GCPnts_AbscissaPoint
460//purpose : rbv for curvilinear parametrization
461//=======================================================================
462
463GCPnts_AbscissaPoint::GCPnts_AbscissaPoint
464 (TheCurve& C,
465 const Standard_Real Abscissa,
466 const Standard_Real U0,
467 const Standard_Real Ui,
468 const Standard_Real Tol)
469{
470 Standard_Real Abscis = Abscissa;
471 Standard_Real UU0 = U0;
472 Standard_Real UUi = Ui;
473 AdvCompute(myComputer, C, Abscis, UU0, UUi, Tol);
474}