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