0027117: BRepClass3d_SolidClassifier doesn't take into account vertex/edge/face toler...
[occt.git] / src / ProjLib / ProjLib_CompProjectedCurve.cxx
CommitLineData
b311480e 1// Created on: 1997-09-23
2// Created by: Roman BORISOV
3// Copyright (c) 1997-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 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
973c2be1 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
42cf5bc1 17
18#include <Adaptor2d_HCurve2d.hxx>
19#include <Adaptor3d_HCurve.hxx>
20#include <Adaptor3d_HSurface.hxx>
7fd59977 21#include <Extrema_ExtCS.hxx>
42cf5bc1 22#include <Extrema_ExtPS.hxx>
7fd59977 23#include <Extrema_GenLocateExtPS.hxx>
7fd59977 24#include <Extrema_POnCurv.hxx>
42cf5bc1 25#include <Extrema_POnSurf.hxx>
7fd59977 26#include <GeomAbs_CurveType.hxx>
27#include <GeomLib.hxx>
42cf5bc1 28#include <gp_Mat2d.hxx>
29#include <gp_Pnt2d.hxx>
30#include <gp_Vec2d.hxx>
31#include <gp_XY.hxx>
32#include <Precision.hxx>
33#include <ProjLib_CompProjectedCurve.hxx>
34#include <ProjLib_HCompProjectedCurve.hxx>
35#include <ProjLib_PrjResolve.hxx>
36#include <Standard_DomainError.hxx>
37#include <Standard_NoSuchObject.hxx>
38#include <Standard_NotImplemented.hxx>
39#include <Standard_OutOfRange.hxx>
40#include <TColgp_HSequenceOfPnt.hxx>
7fd59977 41
7fd59977 42#define FuncTol 1.e-10
43
0797d9d3 44#ifdef OCCT_DEBUG_CHRONO
7fd59977 45#include <OSD_Timer.hxx>
46
47static OSD_Chronometer chr_init_point, chr_dicho_bound;
48
49Standard_EXPORT Standard_Real t_init_point, t_dicho_bound;
50Standard_EXPORT Standard_Integer init_point_count, dicho_bound_count;
51
52static void InitChron(OSD_Chronometer& ch)
53{
6e0fd076 54 ch.Reset();
55 ch.Start();
7fd59977 56}
57
58static void ResultChron( OSD_Chronometer & ch, Standard_Real & time)
59{
6e0fd076 60 Standard_Real tch ;
61 ch.Stop();
62 ch.Show(tch);
63 time=time +tch;
7fd59977 64}
65#endif
66
7fd59977 67
68//=======================================================================
69//function : d1
70//purpose : computes first derivative of the projected curve
71//=======================================================================
72
73static void d1(const Standard_Real t,
6e0fd076 74 const Standard_Real u,
75 const Standard_Real v,
76 gp_Vec2d& V,
77 const Handle(Adaptor3d_HCurve)& Curve,
78 const Handle(Adaptor3d_HSurface)& Surface)
7fd59977 79{
80 gp_Pnt S, C;
81 gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v, DC1_t;
82 Surface->D2(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv);
83 Curve->D1(t, C, DC1_t);
84 gp_Vec Ort(C, S);// Ort = S - C
85
86 gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
87 gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u,
6e0fd076 88 DS1_u*DS1_v + Ort*DS2_uv);
7fd59977 89 gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv,
6e0fd076 90 DS1_v*DS1_v + Ort*DS2_v);
7fd59977 91
92 Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
93 if (fabs(det) < gp::Resolution()) Standard_ConstructionError::Raise();
6e0fd076 94
7fd59977 95 gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det),
6e0fd076 96 gp_XY(-dE_dv.X()/det, dE_du.X()/det));
7fd59977 97
98 V = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
99}
100
101//=======================================================================
102//function : d2
103//purpose : computes second derivative of the projected curve
104//=======================================================================
105
6e0fd076 106static void d2(const Standard_Real t,
107 const Standard_Real u,
108 const Standard_Real v,
109 gp_Vec2d& V1, gp_Vec2d& V2,
110 const Handle(Adaptor3d_HCurve)& Curve,
111 const Handle(Adaptor3d_HSurface)& Surface)
7fd59977 112{
113 gp_Pnt S, C;
114 gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v,
6e0fd076 115 DS3_u, DS3_v, DS3_uuv, DS3_uvv,
116 DC1_t, DC2_t;
7fd59977 117 Surface->D3(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv,
6e0fd076 118 DS3_u, DS3_v, DS3_uuv, DS3_uvv);
7fd59977 119 Curve->D2(t, C, DC1_t, DC2_t);
120 gp_Vec Ort(C, S);
121
122 gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
123 gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u,
6e0fd076 124 DS1_u*DS1_v + Ort*DS2_uv);
7fd59977 125 gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv,
6e0fd076 126 DS1_v*DS1_v + Ort*DS2_v);
7fd59977 127
128 Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
129 if (fabs(det) < gp::Resolution()) Standard_ConstructionError::Raise();
130
131 gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det),
6e0fd076 132 gp_XY(-dE_dv.X()/det, dE_du.X()/det));
7fd59977 133
134 // First derivative
135 V1 = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
136
137 /* Second derivative */
138
139 // Computation of d2E_dt2 = S1
140 gp_Vec2d d2E_dt(-DC2_t*DS1_u, -DC2_t*DS1_v);
141
142 // Computation of 2*(d2E/dtdX)(dX/dt) = S2
143 gp_Vec2d d2E1_dtdX(-DC1_t*DS2_u,
6e0fd076 144 -DC1_t*DS2_uv);
7fd59977 145 gp_Vec2d d2E2_dtdX(-DC1_t*DS2_uv,
6e0fd076 146 -DC1_t*DS2_v);
7fd59977 147 gp_Vec2d S2 = 2*gp_Vec2d(d2E1_dtdX*V1, d2E2_dtdX*V1);
148
149 // Computation of (d2E/dX2)*(dX/dt)2 = S3
150
151 // Row11 = (d2E1/du2, d2E1/dudv)
152 Standard_Real tmp;
153 gp_Vec2d Row11(3*DS1_u*DS2_u + Ort*DS3_u,
6e0fd076 154 tmp = 2*DS1_u*DS2_uv +
155 DS1_v*DS2_u + Ort*DS3_uuv);
7fd59977 156
157 // Row12 = (d2E1/dudv, d2E1/dv2)
158 gp_Vec2d Row12(tmp, DS2_v*DS1_u + 2*DS1_v*DS2_uv +
6e0fd076 159 Ort*DS3_uvv);
7fd59977 160
161 // Row21 = (d2E2/du2, d2E2/dudv)
162 gp_Vec2d Row21(DS2_u*DS1_v + 2*DS1_u*DS2_uv + Ort*DS3_uuv,
6e0fd076 163 tmp = 2*DS2_uv*DS1_v + DS1_u*DS2_v + Ort*DS3_uvv);
7fd59977 164
165 // Row22 = (d2E2/duv, d2E2/dvdv)
166 gp_Vec2d Row22(tmp, 3*DS1_v*DS2_v + Ort*DS3_v);
167
168 gp_Vec2d S3(V1*gp_Vec2d(Row11*V1, Row12*V1),
6e0fd076 169 V1*gp_Vec2d(Row21*V1, Row22*V1));
7fd59977 170
171 gp_Vec2d Sum = d2E_dt + S2 + S3;
172
173 V2 = - gp_Vec2d(gp_Vec2d(M.Row(1))*Sum, gp_Vec2d(M.Row(2))*Sum);
174}
175//=======================================================================
176//function : d1CurveOnSurf
177//purpose : computes first derivative of the 3d projected curve
178//=======================================================================
179
41194117 180#if 0
7fd59977 181static void d1CurvOnSurf(const Standard_Real t,
6e0fd076 182 const Standard_Real u,
183 const Standard_Real v,
184 gp_Vec& V,
185 const Handle(Adaptor3d_HCurve)& Curve,
186 const Handle(Adaptor3d_HSurface)& Surface)
7fd59977 187{
188 gp_Pnt S, C;
189 gp_Vec2d V2d;
190 gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v, DC1_t;
191 Surface->D2(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv);
192 Curve->D1(t, C, DC1_t);
193 gp_Vec Ort(C, S);// Ort = S - C
194
195 gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
196 gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u,
6e0fd076 197 DS1_u*DS1_v + Ort*DS2_uv);
7fd59977 198 gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv,
6e0fd076 199 DS1_v*DS1_v + Ort*DS2_v);
7fd59977 200
201 Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
202 if (fabs(det) < gp::Resolution()) Standard_ConstructionError::Raise();
6e0fd076 203
7fd59977 204 gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det),
6e0fd076 205 gp_XY(-dE_dv.X()/det, dE_du.X()/det));
7fd59977 206
207 V2d = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
208
209 V = DS1_u * V2d.X() + DS1_v * V2d.Y();
210
211}
212#endif
213
214//=======================================================================
215//function : d2CurveOnSurf
216//purpose : computes second derivative of the 3D projected curve
217//=======================================================================
218
6e0fd076 219static void d2CurvOnSurf(const Standard_Real t,
220 const Standard_Real u,
221 const Standard_Real v,
222 gp_Vec& V1 , gp_Vec& V2 ,
223 const Handle(Adaptor3d_HCurve)& Curve,
224 const Handle(Adaptor3d_HSurface)& Surface)
7fd59977 225{
226 gp_Pnt S, C;
227 gp_Vec2d V12d,V22d;
228 gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v,
6e0fd076 229 DS3_u, DS3_v, DS3_uuv, DS3_uvv,
230 DC1_t, DC2_t;
7fd59977 231 Surface->D3(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv,
6e0fd076 232 DS3_u, DS3_v, DS3_uuv, DS3_uvv);
7fd59977 233 Curve->D2(t, C, DC1_t, DC2_t);
234 gp_Vec Ort(C, S);
235
236 gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
237 gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u,
6e0fd076 238 DS1_u*DS1_v + Ort*DS2_uv);
7fd59977 239 gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv,
6e0fd076 240 DS1_v*DS1_v + Ort*DS2_v);
7fd59977 241
242 Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
243 if (fabs(det) < gp::Resolution()) Standard_ConstructionError::Raise();
244
245 gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det),
6e0fd076 246 gp_XY(-dE_dv.X()/det, dE_du.X()/det));
7fd59977 247
248 // First derivative
249 V12d = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
250
251 /* Second derivative */
252
253 // Computation of d2E_dt2 = S1
254 gp_Vec2d d2E_dt(-DC2_t*DS1_u, -DC2_t*DS1_v);
255
256 // Computation of 2*(d2E/dtdX)(dX/dt) = S2
257 gp_Vec2d d2E1_dtdX(-DC1_t*DS2_u,
6e0fd076 258 -DC1_t*DS2_uv);
7fd59977 259 gp_Vec2d d2E2_dtdX(-DC1_t*DS2_uv,
6e0fd076 260 -DC1_t*DS2_v);
7fd59977 261 gp_Vec2d S2 = 2*gp_Vec2d(d2E1_dtdX*V12d, d2E2_dtdX*V12d);
262
263 // Computation of (d2E/dX2)*(dX/dt)2 = S3
264
265 // Row11 = (d2E1/du2, d2E1/dudv)
266 Standard_Real tmp;
267 gp_Vec2d Row11(3*DS1_u*DS2_u + Ort*DS3_u,
6e0fd076 268 tmp = 2*DS1_u*DS2_uv +
269 DS1_v*DS2_u + Ort*DS3_uuv);
7fd59977 270
271 // Row12 = (d2E1/dudv, d2E1/dv2)
272 gp_Vec2d Row12(tmp, DS2_v*DS1_u + 2*DS1_v*DS2_uv +
6e0fd076 273 Ort*DS3_uvv);
7fd59977 274
275 // Row21 = (d2E2/du2, d2E2/dudv)
276 gp_Vec2d Row21(DS2_u*DS1_v + 2*DS1_u*DS2_uv + Ort*DS3_uuv,
6e0fd076 277 tmp = 2*DS2_uv*DS1_v + DS1_u*DS2_v + Ort*DS3_uvv);
7fd59977 278
279 // Row22 = (d2E2/duv, d2E2/dvdv)
280 gp_Vec2d Row22(tmp, 3*DS1_v*DS2_v + Ort*DS3_v);
281
282 gp_Vec2d S3(V12d*gp_Vec2d(Row11*V12d, Row12*V12d),
6e0fd076 283 V12d*gp_Vec2d(Row21*V12d, Row22*V12d));
7fd59977 284
285 gp_Vec2d Sum = d2E_dt + S2 + S3;
286
287 V22d = - gp_Vec2d(gp_Vec2d(M.Row(1))*Sum, gp_Vec2d(M.Row(2))*Sum);
288
289 V1 = DS1_u * V12d.X() + DS1_v * V12d.Y();
290 V2 = DS2_u * V12d.X() *V12d.X()
6e0fd076 291 + DS1_u * V22d.X()
292 + 2 * DS2_uv * V12d.X() *V12d.Y()
293 + DS2_v * V12d.Y() * V12d.Y()
294 + DS1_v * V22d.Y();
7fd59977 295}
296
297//=======================================================================
298//function : ExactBound
299//purpose : computes exact boundary point
300//=======================================================================
301
302static Standard_Boolean ExactBound(gp_Pnt& Sol,
6e0fd076 303 const Standard_Real NotSol,
304 const Standard_Real Tol,
305 const Standard_Real TolU,
306 const Standard_Real TolV,
307 const Handle(Adaptor3d_HCurve)& Curve,
308 const Handle(Adaptor3d_HSurface)& Surface)
7fd59977 309{
310 Standard_Real U0, V0, t, t1, t2, FirstU, LastU, FirstV, LastV;
311 gp_Pnt2d POnS;
312 U0 = Sol.Y();
313 V0 = Sol.Z();
314 FirstU = Surface->FirstUParameter();
315 LastU = Surface->LastUParameter();
316 FirstV = Surface->FirstVParameter();
317 LastV = Surface->LastVParameter();
318 // Here we have to compute the boundary that projection is going to intersect
319 gp_Vec2d D2d;
320 //these variables are to estimate which boundary has more apportunity
321 //to be intersected
322 Standard_Real RU1, RU2, RV1, RV2;
323 d1(Sol.X(), U0, V0, D2d, Curve, Surface);
324 // Here we assume that D2d != (0, 0)
325 if(Abs(D2d.X()) < gp::Resolution())
326 {
327 RU1 = Precision::Infinite();
328 RU2 = Precision::Infinite();
329 RV1 = V0 - FirstV;
330 RV2 = LastV - V0;
331 }
332 else if(Abs(D2d.Y()) < gp::Resolution())
333 {
334 RU1 = U0 - FirstU;
335 RU2 = LastU - U0;
336 RV1 = Precision::Infinite();
337 RV2 = Precision::Infinite();
338 }
339 else
340 {
341 RU1 = gp_Pnt2d(U0, V0).
6e0fd076 342 Distance(gp_Pnt2d(FirstU, V0 + (FirstU - U0)*D2d.Y()/D2d.X()));
7fd59977 343 RU2 = gp_Pnt2d(U0, V0).
6e0fd076 344 Distance(gp_Pnt2d(LastU, V0 + (LastU - U0)*D2d.Y()/D2d.X()));
7fd59977 345 RV1 = gp_Pnt2d(U0, V0).
6e0fd076 346 Distance(gp_Pnt2d(U0 + (FirstV - V0)*D2d.X()/D2d.Y(), FirstV));
7fd59977 347 RV2 = gp_Pnt2d(U0, V0).
6e0fd076 348 Distance(gp_Pnt2d(U0 + (LastV - V0)*D2d.X()/D2d.Y(), LastV));
7fd59977 349 }
350 TColgp_SequenceOfPnt Seq;
351 Seq.Append(gp_Pnt(FirstU, RU1, 2));
352 Seq.Append(gp_Pnt(LastU, RU2, 2));
353 Seq.Append(gp_Pnt(FirstV, RV1, 3));
354 Seq.Append(gp_Pnt(LastV, RV2, 3));
355 Standard_Integer i, j;
356 for(i = 1; i <= 3; i++)
357 for(j = 1; j <= 4-i; j++)
358 if(Seq(j).Y() < Seq(j+1).Y())
359 {
6e0fd076 360 gp_Pnt swp;
361 swp = Seq.Value(j+1);
362 Seq.ChangeValue(j+1) = Seq.Value(j);
363 Seq.ChangeValue(j) = swp;
7fd59977 364 }
365
6e0fd076 366 t = Sol.X();
367 t1 = Min(Sol.X(), NotSol);
368 t2 = Max(Sol.X(), NotSol);
7fd59977 369
6e0fd076 370 Standard_Boolean isDone = Standard_False;
371 while (!Seq.IsEmpty())
372 {
373 gp_Pnt P;
374 P = Seq.Last();
375 Seq.Remove(Seq.Length());
376 ProjLib_PrjResolve aPrjPS(Curve->Curve(),
377 Surface->Surface(),
378 Standard_Integer(P.Z()));
379 if(Standard_Integer(P.Z()) == 2)
380 {
381 aPrjPS.Perform(t, P.X(), V0, gp_Pnt2d(Tol, TolV),
382 gp_Pnt2d(t1, Surface->FirstVParameter()),
383 gp_Pnt2d(t2, Surface->LastVParameter()), FuncTol);
384 if(!aPrjPS.IsDone()) continue;
385 POnS = aPrjPS.Solution();
386 Sol = gp_Pnt(POnS.X(), P.X(), POnS.Y());
387 isDone = Standard_True;
388 break;
389 }
390 else
391 {
392 aPrjPS.Perform(t, U0, P.X(), gp_Pnt2d(Tol, TolU),
393 gp_Pnt2d(t1, Surface->FirstUParameter()),
394 gp_Pnt2d(t2, Surface->LastUParameter()), FuncTol);
395 if(!aPrjPS.IsDone()) continue;
396 POnS = aPrjPS.Solution();
397 Sol = gp_Pnt(POnS.X(), POnS.Y(), P.X());
398 isDone = Standard_True;
399 break;
400 }
401 }
7fd59977 402
6e0fd076 403 return isDone;
7fd59977 404}
405
406//=======================================================================
407//function : DichExactBound
408//purpose : computes exact boundary point
409//=======================================================================
410
411static void DichExactBound(gp_Pnt& Sol,
6e0fd076 412 const Standard_Real NotSol,
413 const Standard_Real Tol,
414 const Standard_Real TolU,
415 const Standard_Real TolV,
416 const Handle(Adaptor3d_HCurve)& Curve,
417 const Handle(Adaptor3d_HSurface)& Surface)
7fd59977 418{
0797d9d3 419#ifdef OCCT_DEBUG_CHRONO
7fd59977 420 InitChron(chr_dicho_bound);
421#endif
422
423 Standard_Real U0, V0, t;
424 gp_Pnt2d POnS;
425 U0 = Sol.Y();
426 V0 = Sol.Z();
427 ProjLib_PrjResolve aPrjPS(Curve->Curve(), Surface->Surface(), 1);
428
429 Standard_Real aNotSol = NotSol;
430 while (fabs(Sol.X() - aNotSol) > Tol)
431 {
432 t = (Sol.X() + aNotSol)/2;
433 aPrjPS.Perform(t, U0, V0, gp_Pnt2d(TolU, TolV),
6e0fd076 434 gp_Pnt2d(Surface->FirstUParameter(),Surface->FirstVParameter()),
435 gp_Pnt2d(Surface->LastUParameter(),Surface->LastVParameter()),
436 FuncTol, Standard_True);
7fd59977 437
438 if (aPrjPS.IsDone())
439 {
440 POnS = aPrjPS.Solution();
441 Sol = gp_Pnt(t, POnS.X(), POnS.Y());
442 U0=Sol.Y();
443 V0=Sol.Z();
444 }
445 else aNotSol = t;
446 }
0797d9d3 447#ifdef OCCT_DEBUG_CHRONO
6e0fd076 448 ResultChron(chr_dicho_bound,t_dicho_bound);
449 dicho_bound_count++;
7fd59977 450#endif
451}
452
453//=======================================================================
454//function : InitialPoint
455//purpose :
456//=======================================================================
457
458static Standard_Boolean InitialPoint(const gp_Pnt& Point,
6e0fd076 459 const Standard_Real t,
460 const Handle(Adaptor3d_HCurve)& C,
461 const Handle(Adaptor3d_HSurface)& S,
462 const Standard_Real TolU,
463 const Standard_Real TolV,
464 Standard_Real& U,
465 Standard_Real& V)
7fd59977 466{
467
6e0fd076 468 ProjLib_PrjResolve aPrjPS(C->Curve(), S->Surface(), 1);
469 Standard_Real ParU,ParV;
470 Extrema_ExtPS aExtPS;
471 aExtPS.Initialize(S->Surface(), S->FirstUParameter(),
472 S->LastUParameter(), S->FirstVParameter(),
473 S->LastVParameter(), TolU, TolV);
7fd59977 474
6e0fd076 475 aExtPS.Perform(Point);
476 Standard_Integer argmin = 0;
477 if (aExtPS.IsDone() && aExtPS.NbExt())
478 {
479 Standard_Integer i, Nend;
480 // Search for the nearest solution which is also a normal projection
481 Nend = aExtPS.NbExt();
482 for(i = 1; i <= Nend; i++)
7fd59977 483 {
6e0fd076 484 Extrema_POnSurf POnS = aExtPS.Point(i);
485 POnS.Parameter(ParU, ParV);
486 aPrjPS.Perform(t, ParU, ParV, gp_Pnt2d(TolU, TolV),
487 gp_Pnt2d(S->FirstUParameter(), S->FirstVParameter()),
488 gp_Pnt2d(S->LastUParameter(), S->LastVParameter()),
489 FuncTol, Standard_True);
490 if(aPrjPS.IsDone() )
491 if (argmin == 0 || aExtPS.SquareDistance(i) < aExtPS.SquareDistance(argmin)) argmin = i;
7fd59977 492 }
6e0fd076 493 }
494 if( argmin == 0 ) return Standard_False;
495 else
496 {
497 Extrema_POnSurf POnS = aExtPS.Point(argmin);
498 POnS.Parameter(U, V);
499 return Standard_True;
500 }
7fd59977 501}
502
503//=======================================================================
504//function : ProjLib_CompProjectedCurve
505//purpose :
506//=======================================================================
507
6e0fd076 508ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve()
cbff1e55 509: myNbCurves(0),
510 myTolU (0.0),
511 myTolV (0.0),
512 myMaxDist (0.0)
7fd59977 513{
514}
515
516//=======================================================================
517//function : ProjLib_CompProjectedCurve
518//purpose :
519//=======================================================================
520
cbff1e55 521ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
522 (const Handle(Adaptor3d_HSurface)& theSurface,
523 const Handle(Adaptor3d_HCurve)& theCurve,
524 const Standard_Real theTolU,
525 const Standard_Real theTolV)
526: mySurface (theSurface),
527 myCurve (theCurve),
528 myNbCurves(0),
529 mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
530 myTolU (theTolU),
531 myTolV (theTolV),
532 myMaxDist (-1.0)
7fd59977 533{
7fd59977 534 Init();
535}
536
537//=======================================================================
538//function : ProjLib_CompProjectedCurve
539//purpose :
540//=======================================================================
541
cbff1e55 542ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
543 (const Handle(Adaptor3d_HSurface)& theSurface,
544 const Handle(Adaptor3d_HCurve)& theCurve,
545 const Standard_Real theTolU,
546 const Standard_Real theTolV,
547 const Standard_Real theMaxDist)
548: mySurface (theSurface),
549 myCurve (theCurve),
550 myNbCurves(0),
551 mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
552 myTolU (theTolU),
553 myTolV (theTolV),
554 myMaxDist (theMaxDist)
7fd59977 555{
7fd59977 556 Init();
557}
558
559//=======================================================================
560//function : Init
561//purpose :
562//=======================================================================
563
6e0fd076 564void ProjLib_CompProjectedCurve::Init()
7fd59977 565{
41194117 566 myTabInt.Nullify();
7fd59977 567
568 Standard_Real Tol;// Tolerance for ExactBound
569 Standard_Integer i, Nend = 0;
570 Standard_Boolean FromLastU=Standard_False;
571
572 //new part (to discard far solutions)
7fd59977 573 Standard_Real TolC = Precision::Confusion(), TolS = Precision::Confusion();
574 Extrema_ExtCS CExt(myCurve->Curve(),
6e0fd076 575 mySurface->Surface(),
576 TolC,
577 TolS);
7fd59977 578 if (CExt.IsDone() && CExt.NbExt())
579 {
6e0fd076 580 // Search for the minimum solution
581 Nend = CExt.NbExt();
aa9d6bec 582 if(myMaxDist > 0 &&
583 // Avoid usage of extrema result that can be wrong for extrusion
584 mySurface->GetType() != GeomAbs_SurfaceOfExtrusion)
6e0fd076 585 {
586 Standard_Real min_val2;
587 min_val2 = CExt.SquareDistance(1);
588 for(i = 2; i <= Nend; i++)
aa9d6bec 589 if (CExt.SquareDistance(i) < min_val2) min_val2 = CExt.SquareDistance(i);
590 if (min_val2 > myMaxDist * myMaxDist)
591 return;
6e0fd076 592 }
593 }
594 // end of new part
7fd59977 595
d1db9125 596 Standard_Real FirstU, LastU, Step, SearchStep, WalkStep, t;
6e0fd076 597
7fd59977 598 FirstU = myCurve->FirstParameter();
599 LastU = myCurve->LastParameter();
d1db9125 600 const Standard_Real GlobalMinStep = 1.e-4;
601 //<GlobalMinStep> is sufficiently small to provide solving from initial point
602 //and, on the other hand, it is sufficiently large to avoid too close solutions.
7fd59977 603 const Standard_Real MinStep = 0.01*(LastU - FirstU),
6e0fd076 604 MaxStep = 0.1*(LastU - FirstU);
7fd59977 605 SearchStep = 10*MinStep;
606 Step = SearchStep;
6e0fd076 607
7fd59977 608 //Initialization of aPrjPS
609 Standard_Real Uinf = mySurface->FirstUParameter();
610 Standard_Real Usup = mySurface->LastUParameter();
611 Standard_Real Vinf = mySurface->FirstVParameter();
612 Standard_Real Vsup = mySurface->LastVParameter();
613
614 ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
615
616 t = FirstU;
617 Standard_Boolean new_part;
618 Standard_Real prevDeb=0.;
619 Standard_Boolean SameDeb=Standard_False;
6e0fd076 620
621
7fd59977 622 gp_Pnt Triple, prevTriple;
623
0d1536ad 624 //Basic loop
7fd59977 625 while(t <= LastU)
626 {
db2a696d 627 // Search for the beginning of a new continuous part
628 // to avoid infinite computation in some difficult cases.
7fd59977 629 new_part = Standard_False;
630 if(t > FirstU && Abs(t-prevDeb) <= Precision::PConfusion()) SameDeb=Standard_True;
631 while(t <= LastU && !new_part && !FromLastU && !SameDeb)
632 {
633 prevDeb=t;
634 if (t == LastU) FromLastU=Standard_True;
635 Standard_Boolean initpoint=Standard_False;
1d47d8d0 636 Standard_Real U = 0., V = 0.;
7fd59977 637 gp_Pnt CPoint;
638 Standard_Real ParT,ParU,ParV;
639
db2a696d 640 // Search an initial point in the list of Extrema Curve-Surface
7fd59977 641 if(Nend != 0 && !CExt.IsParallel())
642 {
6e0fd076 643 for (i=1;i<=Nend;i++)
644 {
645 Extrema_POnCurv P1;
646 Extrema_POnSurf P2;
647 CExt.Points(i,P1,P2);
648 ParT=P1.Parameter();
649 P2.Parameter(ParU, ParV);
650
651 aPrjPS.Perform(ParT, ParU, ParV, gp_Pnt2d(myTolU, myTolV),
652 gp_Pnt2d(mySurface->FirstUParameter(),mySurface->FirstVParameter()),
653 gp_Pnt2d(mySurface->LastUParameter(), mySurface->LastVParameter()),
654 FuncTol, Standard_True);
655 if ( aPrjPS.IsDone() && P1.Parameter() > Max(FirstU,t-Step+Precision::PConfusion())
656 && P1.Parameter() <= t)
657 {
658 t=ParT;
659 U=ParU;
660 V=ParV;
661 CPoint=P1.Value();
662 initpoint = Standard_True;
663 break;
664 }
665 }
7fd59977 666 }
667 if (!initpoint)
668 {
6e0fd076 669 myCurve->D0(t,CPoint);
0797d9d3 670#ifdef OCCT_DEBUG_CHRONO
6e0fd076 671 InitChron(chr_init_point);
7fd59977 672#endif
0d1536ad 673 // PConfusion - use geometric tolerances in extrema / optimization.
674 initpoint=InitialPoint(CPoint, t,myCurve,mySurface, Precision::PConfusion(), Precision::PConfusion(), U, V);
0797d9d3 675#ifdef OCCT_DEBUG_CHRONO
6e0fd076 676 ResultChron(chr_init_point,t_init_point);
677 init_point_count++;
7fd59977 678#endif
6e0fd076 679 }
7fd59977 680 if(initpoint)
681 {
682 // When U or V lie on surface joint in some cases we cannot use them
683 // as initial point for aPrjPS, so we switch them
6e0fd076 684 gp_Vec2d D;
685
d1db9125 686 if ((mySurface->IsUPeriodic() &&
687 Abs(Usup - Uinf - mySurface->UPeriod()) < Precision::Confusion()) ||
688 (mySurface->IsVPeriodic() &&
689 Abs(Vsup - Vinf - mySurface->VPeriod()) < Precision::Confusion()))
6e0fd076 690 {
d1db9125 691 if((Abs(U - Uinf) < mySurface->UResolution(Precision::PConfusion())) &&
692 mySurface->IsUPeriodic())
693 {
694 d1(t, U, V, D, myCurve, mySurface);
695 if (D.X() < 0 ) U = Usup;
696 }
697 else if((Abs(U - Usup) < mySurface->UResolution(Precision::PConfusion())) &&
698 mySurface->IsUPeriodic())
699 {
700 d1(t, U, V, D, myCurve, mySurface);
701 if (D.X() > 0) U = Uinf;
702 }
fa6cd915 703
d1db9125 704 if((Abs(V - Vinf) < mySurface->VResolution(Precision::PConfusion())) &&
705 mySurface->IsVPeriodic())
706 {
707 d1(t, U, V, D, myCurve, mySurface);
708 if (D.Y() < 0) V = Vsup;
709 }
710 else if((Abs(V - Vsup) <= mySurface->VResolution(Precision::PConfusion())) &&
711 mySurface->IsVPeriodic())
712 {
713 d1(t, U, V, D, myCurve, mySurface);
714 if (D.Y() > 0) V = Vinf;
715 }
6e0fd076 716 }
7fd59977 717
718
6e0fd076 719 if (myMaxDist > 0)
7fd59977 720 {
721 // Here we are going to stop if the distance between projection and
722 // corresponding curve point is greater than myMaxDist
6e0fd076 723 gp_Pnt POnS;
724 Standard_Real d;
725 mySurface->D0(U, V, POnS);
726 d = CPoint.Distance(POnS);
727 if (d > myMaxDist)
7fd59977 728 {
6e0fd076 729 mySequence->Clear();
730 myNbCurves = 0;
731 return;
732 }
7fd59977 733 }
6e0fd076 734 Triple = gp_Pnt(t, U, V);
735 if (t != FirstU)
7fd59977 736 {
6e0fd076 737 //Search for exact boundary point
738 Tol = Min(myTolU, myTolV);
51740958 739 gp_Vec2d aD;
740 d1(Triple.X(), Triple.Y(), Triple.Z(), aD, myCurve, mySurface);
741 Tol /= Max(Abs(aD.X()), Abs(aD.Y()));
6e0fd076 742
743 if(!ExactBound(Triple, t - Step, Tol,
744 myTolU, myTolV, myCurve, mySurface))
7fd59977 745 {
0797d9d3 746#ifdef OCCT_DEBUG
6e0fd076 747 cout<<"There is a problem with ExactBound computation"<<endl;
7fd59977 748#endif
6e0fd076 749 DichExactBound(Triple, t - Step, Tol, myTolU, myTolV,
750 myCurve, mySurface);
751 }
752 }
753 new_part = Standard_True;
7fd59977 754 }
755 else
756 {
757 if(t == LastU) break;
758 t += Step;
6e0fd076 759 if(t>LastU)
760 {
761 Step =Step+LastU-t;
762 t=LastU;
763 }
7fd59977 764 }
765 }
766 if (!new_part) break;
767
768
769 //We have found a new continuous part
770 Handle(TColgp_HSequenceOfPnt) hSeq = new TColgp_HSequenceOfPnt();
771 mySequence->Append(hSeq);
772 myNbCurves++;
773 mySequence->Value(myNbCurves)->Append(Triple);
774 prevTriple = Triple;
775
776 if (Triple.X() == LastU) break;//return;
777
778 //Computation of WalkStep
779 gp_Vec D1, D2;
780 Standard_Real MagnD1, MagnD2;
781 d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
782 MagnD1 = D1.Magnitude();
783 MagnD2 = D2.Magnitude();
784 if(MagnD2 < Precision::Confusion()) WalkStep = MaxStep;
785 else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
6e0fd076 786
7fd59977 787 Step = WalkStep;
7fd59977 788
789 t = Triple.X() + Step;
790 if (t > LastU) t = LastU;
1cdee2a6 791 Standard_Real prevStep = Step;
4f0d73a9 792 Standard_Real U0, V0;
793 gp_Pnt2d aLowBorder(mySurface->FirstUParameter(),mySurface->FirstVParameter());
794 gp_Pnt2d aUppBorder(mySurface->LastUParameter(), mySurface->LastVParameter());
795 gp_Pnt2d aTol(myTolU, myTolV);
7fd59977 796 //Here we are trying to prolong continuous part
797 while (t <= LastU && new_part)
798 {
7fd59977 799
1cdee2a6 800 U0 = Triple.Y() + (Step / prevStep) * (Triple.Y() - prevTriple.Y());
801 V0 = Triple.Z() + (Step / prevStep) * (Triple.Z() - prevTriple.Z());
4f0d73a9 802 // adjust U0 to be in [mySurface->FirstUParameter(),mySurface->LastUParameter()]
803 U0 = Min(Max(U0, aLowBorder.X()), aUppBorder.X());
804 // adjust V0 to be in [mySurface->FirstVParameter(),mySurface->LastVParameter()]
805 V0 = Min(Max(V0, aLowBorder.Y()), aUppBorder.Y());
7fd59977 806
4f0d73a9 807
808 aPrjPS.Perform(t, U0, V0, aTol,
809 aLowBorder, aUppBorder, FuncTol, Standard_True);
7fd59977 810 if(!aPrjPS.IsDone())
811 {
d1db9125 812 if (Step <= GlobalMinStep)
7fd59977 813 {
6e0fd076 814 //Search for exact boundary point
815 Tol = Min(myTolU, myTolV);
816 gp_Vec2d D;
817 d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
818 Tol /= Max(Abs(D.X()), Abs(D.Y()));
819
820 if(!ExactBound(Triple, t, Tol, myTolU, myTolV,
821 myCurve, mySurface))
822 {
0797d9d3 823#ifdef OCCT_DEBUG
6e0fd076 824 cout<<"There is a problem with ExactBound computation"<<endl;
7fd59977 825#endif
6e0fd076 826 DichExactBound(Triple, t, Tol, myTolU, myTolV,
827 myCurve, mySurface);
828 }
829
830 if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
831 mySequence->Value(myNbCurves)->Append(Triple);
832 if((LastU - Triple.X()) < Tol) {t = LastU + 1; break;}//return;
833
834 Step = SearchStep;
835 t = Triple.X() + Step;
836 if (t > (LastU-MinStep/2) )
837 {
838 Step =Step+LastU-t;
839 t = LastU;
840 }
6e0fd076 841 new_part = Standard_False;
842 }
7fd59977 843 else
844 {
6e0fd076 845 // decrease step
d1db9125 846 Standard_Real SaveStep = Step;
847 Step /= 2.;
6e0fd076 848 t = Triple .X() + Step;
849 if (t > (LastU-MinStep/4) )
850 {
851 Step =Step+LastU-t;
d1db9125 852 if (Abs(Step - SaveStep) <= Precision::PConfusion())
853 Step = GlobalMinStep; //to avoid looping
6e0fd076 854 t = LastU;
855 }
7fd59977 856 }
857 }
858 // Go further
859 else
860 {
1cdee2a6 861 prevTriple = Triple;
862 prevStep = Step;
6e0fd076 863 Triple = gp_Pnt(t, aPrjPS.Solution().X(), aPrjPS.Solution().Y());
864
db2a696d 865 // Check for possible local traps.
866 UpdateTripleByTrapCriteria(Triple);
1cdee2a6 867
6e0fd076 868 if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
869 mySequence->Value(myNbCurves)->Append(Triple);
870 if (t == LastU) {t = LastU + 1; break;}//return;
6e0fd076 871 //Computation of WalkStep
872 d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
873 MagnD1 = D1.Magnitude();
874 MagnD2 = D2.Magnitude();
875 if(MagnD2 < Precision::Confusion() ) WalkStep = MaxStep;
876 else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
877
878 Step = WalkStep;
879 t += Step;
880 if (t > (LastU-MinStep/2) )
1cdee2a6 881 {
6e0fd076 882 Step =Step+LastU-t;
883 t = LastU;
884 }
7fd59977 885 }
886 }
887 }
db2a696d 888 // Sequence post-proceeding.
7fd59977 889 Standard_Integer j;
890
6e0fd076 891 // 1. Removing poor parts
7fd59977 892 Standard_Integer NbPart=myNbCurves;
893 Standard_Integer ipart=1;
894 for(i = 1; i <= NbPart; i++) {
6e0fd076 895 // Standard_Integer NbPoints = mySequence->Value(i)->Length();
7fd59977 896 if(mySequence->Value(ipart)->Length() < 2) {
897 mySequence->Remove(ipart);
898 myNbCurves--;
899 }
900 else ipart++;
901 }
902
903 if(myNbCurves == 0) return;
904
6e0fd076 905 // 2. Removing common parts of bounds
7fd59977 906 for(i = 1; i < myNbCurves; i++)
907 {
908 if(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() >=
6e0fd076 909 mySequence->Value(i+1)->Value(1).X())
7fd59977 910 mySequence->ChangeValue(i+1)->ChangeValue(1).SetX(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() + 1.e-12);
911 }
912
6e0fd076 913 // 3. Computation of the maximum distance from each part of curve to surface
7fd59977 914
915 myMaxDistance = new TColStd_HArray1OfReal(1, myNbCurves);
916 myMaxDistance->Init(0);
917 for(i = 1; i <= myNbCurves; i++)
918 for(j = 1; j <= mySequence->Value(i)->Length(); j++)
919 {
51740958 920 gp_Pnt POnC, POnS, aTriple;
7fd59977 921 Standard_Real Distance;
51740958 922 aTriple = mySequence->Value(i)->Value(j);
923 myCurve->D0(aTriple.X(), POnC);
924 mySurface->D0(aTriple.Y(), aTriple.Z(), POnS);
7fd59977 925 Distance = POnC.Distance(POnS);
926 if (myMaxDistance->Value(i) < Distance)
6e0fd076 927 myMaxDistance->ChangeValue(i) = Distance;
7fd59977 928 }
929
930
6e0fd076 931 // 4. Check the projection to be a single point
7fd59977 932
6e0fd076 933 gp_Pnt2d Pmoy, Pcurr, P;
934 Standard_Real AveU, AveV;
935 mySnglPnts = new TColStd_HArray1OfBoolean(1, myNbCurves);
936 for(i = 1; i <= myNbCurves; i++) mySnglPnts->SetValue(i, Standard_True);
7fd59977 937
6e0fd076 938 for(i = 1; i <= myNbCurves; i++)
939 {
940 //compute an average U and V
7fd59977 941
6e0fd076 942 for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++)
943 {
944 AveU += mySequence->Value(i)->Value(j).Y();
945 AveV += mySequence->Value(i)->Value(j).Z();
946 }
947 AveU /= mySequence->Value(i)->Length();
948 AveV /= mySequence->Value(i)->Length();
7fd59977 949
6e0fd076 950 Pmoy.SetCoord(AveU,AveV);
951 for(j = 1; j <= mySequence->Value(i)->Length(); j++)
952 {
953 Pcurr =
954 gp_Pnt2d(mySequence->Value(i)->Value(j).Y(), mySequence->Value(i)->Value(j).Z());
955 if (Pcurr.Distance(Pmoy) > ((myTolU < myTolV) ? myTolV : myTolU))
7fd59977 956 {
6e0fd076 957 mySnglPnts->SetValue(i, Standard_False);
958 break;
959 }
960 }
7fd59977 961 }
7fd59977 962
6e0fd076 963 // 5. Check the projection to be an isoparametric curve of the surface
7fd59977 964
6e0fd076 965 myUIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
966 for(i = 1; i <= myNbCurves; i++) myUIso->SetValue(i, Standard_True);
7fd59977 967
6e0fd076 968 myVIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
969 for(i = 1; i <= myNbCurves; i++) myVIso->SetValue(i, Standard_True);
7fd59977 970
6e0fd076 971 for(i = 1; i <= myNbCurves; i++) {
972 if (IsSinglePnt(i, P)|| mySequence->Value(i)->Length() <=2) {
973 myUIso->SetValue(i, Standard_False);
974 myVIso->SetValue(i, Standard_False);
975 continue;
976 }
7fd59977 977
6e0fd076 978 // new test for isoparametrics
7fd59977 979
6e0fd076 980 if ( mySequence->Value(i)->Length() > 2) {
981 //compute an average U and V
7fd59977 982
6e0fd076 983 for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++) {
984 AveU += mySequence->Value(i)->Value(j).Y();
985 AveV += mySequence->Value(i)->Value(j).Z();
986 }
987 AveU /= mySequence->Value(i)->Length();
988 AveV /= mySequence->Value(i)->Length();
7fd59977 989
6e0fd076 990 // is i-part U-isoparametric ?
991 for(j = 1; j <= mySequence->Value(i)->Length(); j++)
992 {
993 if(Abs(mySequence->Value(i)->Value(j).Y() - AveU) > myTolU)
994 {
995 myUIso->SetValue(i, Standard_False);
996 break;
997 }
998 }
999
1000 // is i-part V-isoparametric ?
1001 for(j = 1; j <= mySequence->Value(i)->Length(); j++)
1002 {
1003 if(Abs(mySequence->Value(i)->Value(j).Z() - AveV) > myTolV)
1004 {
1005 myVIso->SetValue(i, Standard_False);
1006 break;
1007 }
1008 }
1009 //
7fd59977 1010 }
1011 }
7fd59977 1012}
1013//=======================================================================
1014//function : Load
1015//purpose :
1016//=======================================================================
1017
1018void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HSurface)& S)
1019{
1020 mySurface = S;
1021}
1022
1023//=======================================================================
1024//function : Load
1025//purpose :
1026//=======================================================================
1027
1028void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HCurve)& C)
1029{
1030 myCurve = C;
1031}
1032
1033//=======================================================================
1034//function : GetSurface
1035//purpose :
1036//=======================================================================
1037
6e0fd076 1038const Handle(Adaptor3d_HSurface)& ProjLib_CompProjectedCurve::GetSurface() const
7fd59977 1039{
1040 return mySurface;
1041}
1042
1043
1044//=======================================================================
1045//function : GetCurve
1046//purpose :
1047//=======================================================================
1048
6e0fd076 1049const Handle(Adaptor3d_HCurve)& ProjLib_CompProjectedCurve::GetCurve() const
7fd59977 1050{
1051 return myCurve;
1052}
1053
1054//=======================================================================
1055//function : GetTolerance
1056//purpose :
1057//=======================================================================
1058
6e0fd076 1059void ProjLib_CompProjectedCurve::GetTolerance(Standard_Real& TolU,
1060 Standard_Real& TolV) const
7fd59977 1061{
1062 TolU = myTolU;
1063 TolV = myTolV;
1064}
1065
1066//=======================================================================
1067//function : NbCurves
1068//purpose :
1069//=======================================================================
1070
6e0fd076 1071Standard_Integer ProjLib_CompProjectedCurve::NbCurves() const
7fd59977 1072{
1073 return myNbCurves;
1074}
1075//=======================================================================
1076//function : Bounds
1077//purpose :
1078//=======================================================================
1079
6e0fd076 1080void ProjLib_CompProjectedCurve::Bounds(const Standard_Integer Index,
1081 Standard_Real& Udeb,
1082 Standard_Real& Ufin) const
7fd59977 1083{
1084 if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1085 Udeb = mySequence->Value(Index)->Value(1).X();
1086 Ufin = mySequence->Value(Index)->Value(mySequence->Value(Index)->Length()).X();
1087}
1088//=======================================================================
1089//function : IsSinglePnt
1090//purpose :
1091//=======================================================================
1092
6e0fd076 1093Standard_Boolean ProjLib_CompProjectedCurve::IsSinglePnt(const Standard_Integer Index, gp_Pnt2d& P) const
7fd59977 1094{
1095 if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1096 P = gp_Pnt2d(mySequence->Value(Index)->Value(1).Y(), mySequence->Value(Index)->Value(1).Z());
1097 return mySnglPnts->Value(Index);
1098}
1099
1100//=======================================================================
1101//function : IsUIso
1102//purpose :
1103//=======================================================================
1104
6e0fd076 1105Standard_Boolean ProjLib_CompProjectedCurve::IsUIso(const Standard_Integer Index, Standard_Real& U) const
7fd59977 1106{
1107 if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1108 U = mySequence->Value(Index)->Value(1).Y();
1109 return myUIso->Value(Index);
1110}
1111//=======================================================================
1112//function : IsVIso
1113//purpose :
1114//=======================================================================
1115
6e0fd076 1116Standard_Boolean ProjLib_CompProjectedCurve::IsVIso(const Standard_Integer Index, Standard_Real& V) const
7fd59977 1117{
1118 if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1119 V = mySequence->Value(Index)->Value(1).Z();
1120 return myVIso->Value(Index);
1121}
1122//=======================================================================
1123//function : Value
1124//purpose :
1125//=======================================================================
1126
6e0fd076 1127gp_Pnt2d ProjLib_CompProjectedCurve::Value(const Standard_Real t) const
7fd59977 1128{
1129 gp_Pnt2d P;
1130 D0(t, P);
1131 return P;
1132}
1133//=======================================================================
1134//function : D0
1135//purpose :
1136//=======================================================================
1137
6e0fd076 1138void ProjLib_CompProjectedCurve::D0(const Standard_Real U,gp_Pnt2d& P) const
7fd59977 1139{
1140 Standard_Integer i, j;
1141 Standard_Real Udeb, Ufin;
1142 Standard_Boolean found = Standard_False;
1143
1144 for(i = 1; i <= myNbCurves; i++)
1145 {
1146 Bounds(i, Udeb, Ufin);
1147 if (U >= Udeb && U <= Ufin)
1148 {
1149 found = Standard_True;
1150 break;
1151 }
1152 }
1153 if (!found) Standard_DomainError::Raise("ProjLib_CompProjectedCurve::D0");
1154
1155 Standard_Real U0, V0;
1156
1157 Standard_Integer End = mySequence->Value(i)->Length();
1158 for(j = 1; j < End; j++)
1159 if ((U >= mySequence->Value(i)->Value(j).X()) && (U <= mySequence->Value(i)->Value(j + 1).X())) break;
1160
6e0fd076 1161 // U0 = mySequence->Value(i)->Value(j).Y();
1162 // V0 = mySequence->Value(i)->Value(j).Z();
7fd59977 1163
6e0fd076 1164 // Cubic Interpolation
7fd59977 1165 if(mySequence->Value(i)->Length() < 4 ||
1166 (Abs(U-mySequence->Value(i)->Value(j).X()) <= Precision::PConfusion()) )
1167 {
1168 U0 = mySequence->Value(i)->Value(j).Y();
1169 V0 = mySequence->Value(i)->Value(j).Z();
1170 }
1171 else if (Abs(U-mySequence->Value(i)->Value(j+1).X())
6e0fd076 1172 <= Precision::PConfusion())
7fd59977 1173 {
1174 U0 = mySequence->Value(i)->Value(j+1).Y();
1175 V0 = mySequence->Value(i)->Value(j+1).Z();
1176 }
1177 else
1178 {
1179 if (j == 1) j = 2;
1180 if (j > mySequence->Value(i)->Length() - 2)
6e0fd076 1181 j = mySequence->Value(i)->Length() - 2;
1182
7fd59977 1183 gp_Vec2d I1, I2, I3, I21, I22, I31, Y1, Y2, Y3, Y4, Res;
1184 Standard_Real X1, X2, X3, X4;
6e0fd076 1185
7fd59977 1186 X1 = mySequence->Value(i)->Value(j - 1).X();
1187 X2 = mySequence->Value(i)->Value(j).X();
1188 X3 = mySequence->Value(i)->Value(j + 1).X();
1189 X4 = mySequence->Value(i)->Value(j + 2).X();
6e0fd076 1190
7fd59977 1191 Y1 = gp_Vec2d(mySequence->Value(i)->Value(j - 1).Y(),
6e0fd076 1192 mySequence->Value(i)->Value(j - 1).Z());
7fd59977 1193 Y2 = gp_Vec2d(mySequence->Value(i)->Value(j).Y(),
6e0fd076 1194 mySequence->Value(i)->Value(j).Z());
7fd59977 1195 Y3 = gp_Vec2d(mySequence->Value(i)->Value(j + 1).Y(),
6e0fd076 1196 mySequence->Value(i)->Value(j + 1).Z());
7fd59977 1197 Y4 = gp_Vec2d(mySequence->Value(i)->Value(j + 2).Y(),
6e0fd076 1198 mySequence->Value(i)->Value(j + 2).Z());
1199
7fd59977 1200 I1 = (Y1 - Y2)/(X1 - X2);
1201 I2 = (Y2 - Y3)/(X2 - X3);
1202 I3 = (Y3 - Y4)/(X3 - X4);
6e0fd076 1203
7fd59977 1204 I21 = (I1 - I2)/(X1 - X3);
1205 I22 = (I2 - I3)/(X2 - X4);
6e0fd076 1206
7fd59977 1207 I31 = (I21 - I22)/(X1 - X4);
6e0fd076 1208
7fd59977 1209 Res = Y1 + (U - X1)*(I1 + (U - X2)*(I21 + (U - X3)*I31));
6e0fd076 1210
7fd59977 1211 U0 = Res.X();
1212 V0 = Res.Y();
1213
1214 if(U0 < mySurface->FirstUParameter()) U0 = mySurface->FirstUParameter();
1215 else if(U0 > mySurface->LastUParameter()) U0 = mySurface->LastUParameter();
1216
1217 if(V0 < mySurface->FirstVParameter()) V0 = mySurface->FirstVParameter();
1218 else if(V0 > mySurface->LastVParameter()) V0 = mySurface->LastVParameter();
1219 }
1220 //End of cubic interpolation
1221
1222 ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
1223 aPrjPS.Perform(U, U0, V0, gp_Pnt2d(myTolU, myTolV),
6e0fd076 1224 gp_Pnt2d(mySurface->FirstUParameter(), mySurface->FirstVParameter()),
1225 gp_Pnt2d(mySurface->LastUParameter(), mySurface->LastVParameter()));
d1db9125 1226 if (aPrjPS.IsDone())
1227 P = aPrjPS.Solution();
1228 else
1229 {
1230 gp_Pnt thePoint = myCurve->Value(U);
1231 Extrema_ExtPS aExtPS(thePoint, mySurface->Surface(), myTolU, myTolV);
1232 if (aExtPS.IsDone() && aExtPS.NbExt())
1233 {
51740958 1234 Standard_Integer k, Nend, imin = 1;
d1db9125 1235 // Search for the nearest solution which is also a normal projection
1236 Nend = aExtPS.NbExt();
51740958 1237 for(k = 2; k <= Nend; k++)
1238 if (aExtPS.SquareDistance(k) < aExtPS.SquareDistance(imin))
1239 imin = k;
d1db9125 1240 const Extrema_POnSurf& POnS = aExtPS.Point(imin);
1241 Standard_Real ParU,ParV;
1242 POnS.Parameter(ParU, ParV);
1243 P.SetCoord(ParU, ParV);
1244 }
1245 else
1246 P.SetCoord(U0,V0);
1247 }
7fd59977 1248}
1249//=======================================================================
1250//function : D1
1251//purpose :
1252//=======================================================================
1253
6e0fd076 1254void ProjLib_CompProjectedCurve::D1(const Standard_Real t,
1255 gp_Pnt2d& P,
1256 gp_Vec2d& V) const
7fd59977 1257{
1258 Standard_Real u, v;
1259 D0(t, P);
1260 u = P.X();
1261 v = P.Y();
1262 d1(t, u, v, V, myCurve, mySurface);
1263}
1264//=======================================================================
1265//function : D2
1266//purpose :
1267//=======================================================================
1268
6e0fd076 1269void ProjLib_CompProjectedCurve::D2(const Standard_Real t,
1270 gp_Pnt2d& P,
1271 gp_Vec2d& V1,
1272 gp_Vec2d& V2) const
7fd59977 1273{
1274 Standard_Real u, v;
1275 D0(t, P);
1276 u = P.X();
1277 v = P.Y();
1278 d2(t, u, v, V1, V2, myCurve, mySurface);
1279}
1280//=======================================================================
1281//function : DN
1282//purpose :
1283//=======================================================================
1284
1285gp_Vec2d ProjLib_CompProjectedCurve::DN(const Standard_Real t,
6e0fd076 1286 const Standard_Integer N) const
7fd59977 1287{
1288 if (N < 1 ) Standard_OutOfRange::Raise("ProjLib_CompProjectedCurve : N must be greater than 0");
1289 else if (N ==1)
1290 {
6e0fd076 1291 gp_Pnt2d P;
1292 gp_Vec2d V;
1293 D1(t,P,V);
1294 return V;
1295 }
7fd59977 1296 else if ( N==2)
1297 {
6e0fd076 1298 gp_Pnt2d P;
1299 gp_Vec2d V1,V2;
1300 D2(t,P,V1,V2);
1301 return V2;
7fd59977 1302 }
1303 else if (N > 2 )
6e0fd076 1304 Standard_NotImplemented::Raise("ProjLib_CompProjectedCurve::DN");
7fd59977 1305 return gp_Vec2d();
1306}
1307
1308//=======================================================================
1309//function : GetSequence
1310//purpose :
1311//=======================================================================
1312
6e0fd076 1313const Handle(ProjLib_HSequenceOfHSequenceOfPnt)& ProjLib_CompProjectedCurve::GetSequence() const
7fd59977 1314{
1315 return mySequence;
1316}
1317//=======================================================================
1318//function : FirstParameter
1319//purpose :
1320//=======================================================================
1321
6e0fd076 1322Standard_Real ProjLib_CompProjectedCurve::FirstParameter() const
7fd59977 1323{
1324 return myCurve->FirstParameter();
1325}
1326
1327//=======================================================================
1328//function : LastParameter
1329//purpose :
1330//=======================================================================
1331
6e0fd076 1332Standard_Real ProjLib_CompProjectedCurve::LastParameter() const
7fd59977 1333{
1334 return myCurve->LastParameter();
1335}
1336
1337//=======================================================================
1338//function : MaxDistance
1339//purpose :
1340//=======================================================================
1341
6e0fd076 1342Standard_Real ProjLib_CompProjectedCurve::MaxDistance(const Standard_Integer Index) const
7fd59977 1343{
1344 if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1345 return myMaxDistance->Value(Index);
1346}
1347
1348//=======================================================================
1349//function : NbIntervals
1350//purpose :
1351//=======================================================================
1352
6e0fd076 1353Standard_Integer ProjLib_CompProjectedCurve::NbIntervals(const GeomAbs_Shape S) const
7fd59977 1354{
41194117 1355 const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt.Nullify();
7fd59977 1356 BuildIntervals(S);
41194117 1357 return myTabInt->Length() - 1;
7fd59977 1358}
1359
1360//=======================================================================
1361//function : Intervals
1362//purpose :
1363//=======================================================================
1364
6e0fd076 1365void ProjLib_CompProjectedCurve::Intervals(TColStd_Array1OfReal& T,const GeomAbs_Shape S) const
7fd59977 1366{
41194117
K
1367 if (myTabInt.IsNull()) BuildIntervals (S);
1368 T = myTabInt->Array1();
7fd59977 1369}
1370
1371//=======================================================================
1372//function : BuildIntervals
1373//purpose :
1374//=======================================================================
1375
6e0fd076 1376void ProjLib_CompProjectedCurve::BuildIntervals(const GeomAbs_Shape S) const
7fd59977 1377{
7fd59977 1378 GeomAbs_Shape SforS = GeomAbs_CN;
7fd59977 1379 switch(S) {
1380 case GeomAbs_C0:
1381 SforS = GeomAbs_C1;
1382 break;
1383 case GeomAbs_C1:
1384 SforS = GeomAbs_C2;
1385 break;
1386 case GeomAbs_C2:
1387 SforS = GeomAbs_C3;
1388 break;
1389 case GeomAbs_C3:
1390 SforS = GeomAbs_CN;
1391 break;
1392 case GeomAbs_CN:
1393 SforS = GeomAbs_CN;
1394 break;
1395 default:
1396 Standard_OutOfRange::Raise();
1397 }
1398 Standard_Integer i, j, k;
1399 Standard_Integer NbIntCur = myCurve->NbIntervals(S);
1400 Standard_Integer NbIntSurU = mySurface->NbUIntervals(SforS);
1401 Standard_Integer NbIntSurV = mySurface->NbVIntervals(SforS);
1402
1403 TColStd_Array1OfReal CutPntsT(1, NbIntCur+1);
1404 TColStd_Array1OfReal CutPntsU(1, NbIntSurU+1);
1405 TColStd_Array1OfReal CutPntsV(1, NbIntSurV+1);
1406
1407 myCurve->Intervals(CutPntsT, S);
1408 mySurface->UIntervals(CutPntsU, SforS);
1409 mySurface->VIntervals(CutPntsV, SforS);
1410
1411 Standard_Real Tl, Tr, Ul, Ur, Vl, Vr, Tol;
1412
1413 Handle(TColStd_HArray1OfReal) BArr = NULL,
6e0fd076 1414 CArr = NULL,
1415 UArr = NULL,
1416 VArr = NULL;
7fd59977 1417
1418 // proccessing projection bounds
1419 BArr = new TColStd_HArray1OfReal(1, 2*myNbCurves);
1420 for(i = 1; i <= myNbCurves; i++)
1421 Bounds(i, BArr->ChangeValue(2*i - 1), BArr->ChangeValue(2*i));
1422
1423 // proccessing curve discontinuities
1424 if(NbIntCur > 1) {
1425 CArr = new TColStd_HArray1OfReal(1, NbIntCur - 1);
1426 for(i = 1; i <= CArr->Length(); i++)
1427 CArr->ChangeValue(i) = CutPntsT(i + 1);
1428 }
1429
1430 // proccessing U-surface discontinuities
1431 TColStd_SequenceOfReal TUdisc;
1432
1433 for(k = 2; k <= NbIntSurU; k++) {
6e0fd076 1434 // cout<<"CutPntsU("<<k<<") = "<<CutPntsU(k)<<endl;
7fd59977 1435 for(i = 1; i <= myNbCurves; i++)
1436 for(j = 1; j < mySequence->Value(i)->Length(); j++) {
6e0fd076 1437 Ul = mySequence->Value(i)->Value(j).Y();
1438 Ur = mySequence->Value(i)->Value(j + 1).Y();
1439
1440 if(Abs(Ul - CutPntsU(k)) <= myTolU)
1441 TUdisc.Append(mySequence->Value(i)->Value(j).X());
1442 else if(Abs(Ur - CutPntsU(k)) <= myTolU)
1443 TUdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1444 else if((Ul < CutPntsU(k) && CutPntsU(k) < Ur) ||
0ebaa4db 1445 (Ur < CutPntsU(k) && CutPntsU(k) < Ul))
7fd59977 1446 {
6e0fd076 1447 Standard_Real V;
1448 V = (mySequence->Value(i)->Value(j).Z()
7fd59977 1449 + mySequence->Value(i)->Value(j +1).Z())/2;
6e0fd076 1450 ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 2);
1451
1452 gp_Vec2d D;
1453 gp_Pnt Triple;
1454 Triple = mySequence->Value(i)->Value(j);
1455 d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1456 if (Abs(D.X()) < Precision::Confusion())
1457 Tol = myTolU;
1458 else
1459 Tol = Min(myTolU, myTolU / Abs(D.X()));
1460
1461 Tl = mySequence->Value(i)->Value(j).X();
1462 Tr = mySequence->Value(i)->Value(j + 1).X();
1463
1464 Solver.Perform((Tl + Tr)/2, CutPntsU(k), V,
1465 gp_Pnt2d(Tol, myTolV),
1466 gp_Pnt2d(Tl, mySurface->FirstVParameter()),
1467 gp_Pnt2d(Tr, mySurface->LastVParameter()));
1468 //
1469 if(Solver.IsDone())
1470 {
1471 TUdisc.Append(Solver.Solution().X());
1472 }
1473 }
7fd59977 1474 }
1475 }
1476 for(i = 2; i <= TUdisc.Length(); i++)
1477 if(TUdisc(i) - TUdisc(i-1) < Precision::PConfusion())
1478 TUdisc.Remove(i--);
1479
1480 if(TUdisc.Length())
1481 {
1482 UArr = new TColStd_HArray1OfReal(1, TUdisc.Length());
1483 for(i = 1; i <= UArr->Length(); i++)
1484 UArr->ChangeValue(i) = TUdisc(i);
1485 }
1486 // proccessing V-surface discontinuities
1487 TColStd_SequenceOfReal TVdisc;
1488
1489 for(k = 2; k <= NbIntSurV; k++)
1490 for(i = 1; i <= myNbCurves; i++)
1491 {
6e0fd076 1492 // cout<<"CutPntsV("<<k<<") = "<<CutPntsV(k)<<endl;
7fd59977 1493 for(j = 1; j < mySequence->Value(i)->Length(); j++) {
1494
6e0fd076 1495 Vl = mySequence->Value(i)->Value(j).Z();
1496 Vr = mySequence->Value(i)->Value(j + 1).Z();
7fd59977 1497
6e0fd076 1498 if(Abs(Vl - CutPntsV(k)) <= myTolV)
1499 TVdisc.Append(mySequence->Value(i)->Value(j).X());
1500 else if (Abs(Vr - CutPntsV(k)) <= myTolV)
1501 TVdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1502 else if((Vl < CutPntsV(k) && CutPntsV(k) < Vr) ||
0ebaa4db 1503 (Vr < CutPntsV(k) && CutPntsV(k) < Vl))
7fd59977 1504 {
6e0fd076 1505 Standard_Real U;
1506 U = (mySequence->Value(i)->Value(j).Y()
1507 + mySequence->Value(i)->Value(j +1).Y())/2;
1508 ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 3);
1509
1510 gp_Vec2d D;
1511 gp_Pnt Triple;
1512 Triple = mySequence->Value(i)->Value(j);
1513 d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1514 if (Abs(D.Y()) < Precision::Confusion())
1515 Tol = myTolV;
1516 else
1517 Tol = Min(myTolV, myTolV / Abs(D.Y()));
1518
1519 Tl = mySequence->Value(i)->Value(j).X();
1520 Tr = mySequence->Value(i)->Value(j + 1).X();
1521
1522 Solver.Perform((Tl + Tr)/2, U, CutPntsV(k),
1523 gp_Pnt2d(Tol, myTolV),
1524 gp_Pnt2d(Tl, mySurface->FirstUParameter()),
1525 gp_Pnt2d(Tr, mySurface->LastUParameter()));
1526 //
1527 if(Solver.IsDone())
1528 {
1529 TVdisc.Append(Solver.Solution().X());
1530 }
1531 }
7fd59977 1532 }
6e0fd076 1533 }
1534 for(i = 2; i <= TVdisc.Length(); i++)
1535 if(TVdisc(i) - TVdisc(i-1) < Precision::PConfusion())
1536 TVdisc.Remove(i--);
7fd59977 1537
6e0fd076 1538 if(TVdisc.Length())
1539 {
1540 VArr = new TColStd_HArray1OfReal(1, TVdisc.Length());
1541 for(i = 1; i <= VArr->Length(); i++)
1542 VArr->ChangeValue(i) = TVdisc(i);
1543 }
7fd59977 1544
6e0fd076 1545 // fusion
1546 TColStd_SequenceOfReal Fusion;
1547 if(!CArr.IsNull())
1548 {
1549 GeomLib::FuseIntervals(BArr->ChangeArray1(),
1550 CArr->ChangeArray1(),
1551 Fusion, Precision::PConfusion());
1552 BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1553 for(i = 1; i <= BArr->Length(); i++)
1554 BArr->ChangeValue(i) = Fusion(i);
1555 Fusion.Clear();
1556 }
7fd59977 1557
6e0fd076 1558 if(!UArr.IsNull())
1559 {
1560 GeomLib::FuseIntervals(BArr->ChangeArray1(),
1561 UArr->ChangeArray1(),
1562 Fusion, Precision::PConfusion());
1563 BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1564 for(i = 1; i <= BArr->Length(); i++)
1565 BArr->ChangeValue(i) = Fusion(i);
1566 Fusion.Clear();
1567 }
7fd59977 1568
6e0fd076 1569 if(!VArr.IsNull())
1570 {
1571 GeomLib::FuseIntervals(BArr->ChangeArray1(),
1572 VArr->ChangeArray1(),
1573 Fusion, Precision::PConfusion());
1574 BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1575 for(i = 1; i <= BArr->Length(); i++)
1576 BArr->ChangeValue(i) = Fusion(i);
1577 }
7fd59977 1578
6e0fd076 1579 const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt = new TColStd_HArray1OfReal(1, BArr->Length());
1580 for(i = 1; i <= BArr->Length(); i++)
1581 myTabInt->ChangeValue(i) = BArr->Value(i);
7fd59977 1582
1583}
1584
1585//=======================================================================
1586//function : Trim
1587//purpose :
1588//=======================================================================
1589
1590Handle(Adaptor2d_HCurve2d) ProjLib_CompProjectedCurve::Trim
6e0fd076 1591 (const Standard_Real First,
1592 const Standard_Real Last,
1593 const Standard_Real Tol) const
7fd59977 1594{
1595 Handle(ProjLib_HCompProjectedCurve) HCS =
6e0fd076 1596 new ProjLib_HCompProjectedCurve(*this);
7fd59977 1597 HCS->ChangeCurve2d().Load(mySurface);
1598 HCS->ChangeCurve2d().Load(myCurve->Trim(First,Last,Tol));
1599 return HCS;
1600}
1601
1602//=======================================================================
1603//function : GetType
1604//purpose :
1605//=======================================================================
1606
1607GeomAbs_CurveType ProjLib_CompProjectedCurve::GetType() const
1608{
1609 return GeomAbs_OtherCurve;
1610}
db2a696d 1611
1612//=======================================================================
1613//function : UpdateTripleByTrapCriteria
1614//purpose :
1615//=======================================================================
1616void ProjLib_CompProjectedCurve::UpdateTripleByTrapCriteria(gp_Pnt &thePoint) const
1617{
1618 Standard_Boolean isProblemsPossible = Standard_False;
1619 // Check possible traps cases:
1620
1621 // 25892 bug.
1622 if (mySurface->GetType() == GeomAbs_SurfaceOfRevolution)
1623 {
1624 // Compute maximal deviation from 3D and choose the biggest one.
1625 Standard_Real aVRes = mySurface->VResolution(Precision::Confusion());
1626 Standard_Real aMaxTol = Max(Precision::PConfusion(), aVRes);
1627
1628 if (Abs (thePoint.Z() - mySurface->FirstVParameter()) < aMaxTol ||
1629 Abs (thePoint.Z() - mySurface->LastVParameter() ) < aMaxTol )
1630 {
1631 isProblemsPossible = Standard_True;
1632 }
1633 }
1634
1635 // 27135 bug. Trap on degenerated edge.
1636 if (mySurface->GetType() == GeomAbs_Sphere &&
1637 (Abs (thePoint.Z() - mySurface->FirstVParameter()) < Precision::PConfusion() ||
1638 Abs (thePoint.Z() - mySurface->LastVParameter() ) < Precision::PConfusion() ||
1639 Abs (thePoint.Y() - mySurface->FirstUParameter()) < Precision::PConfusion() ||
1640 Abs (thePoint.Y() - mySurface->LastUParameter() ) < Precision::PConfusion() ))
1641 {
1642 isProblemsPossible = Standard_True;
1643 }
1644
1645 if (!isProblemsPossible)
1646 return;
1647
1648 Standard_Real U,V;
0d1536ad 1649 Standard_Boolean isDone =
1650 InitialPoint(myCurve->Value(thePoint.X()), thePoint.X(), myCurve, mySurface,
1651 Precision::PConfusion(), Precision::PConfusion(), U, V);
1652
1653 if (!isDone)
1654 return;
db2a696d 1655
1656 // Restore original position in case of period jump.
1657 if (mySurface->IsUPeriodic() &&
1658 Abs (Abs(U - thePoint.Y()) - mySurface->UPeriod()) < Precision::PConfusion())
1659 {
1660 U = thePoint.Y();
1661 }
1662 if (mySurface->IsVPeriodic() &&
1663 Abs (Abs(V - thePoint.Z()) - mySurface->VPeriod()) < Precision::PConfusion())
1664 {
1665 V = thePoint.Z();
1666 }
1667 thePoint.SetY(U);
1668 thePoint.SetZ(V);
1669}