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