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