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