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