1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
16 #include <Adaptor3d_HSurface.hxx>
17 #include <Adaptor3d_HSurfaceTool.hxx>
18 #include <Extrema_GenLocateExtPS.hxx>
19 #include <Geom_Surface.hxx>
22 #include <gp_Pnt2d.hxx>
23 #include <IntImp_ComputeTangence.hxx>
24 #include <IntSurf_LineOn2S.hxx>
25 #include <IntSurf_PntOn2S.hxx>
26 #include <IntWalk_PWalking.hxx>
27 #include <IntWalk_StatusDeflection.hxx>
28 #include <math_FunctionSetRoot.hxx>
29 #include <Precision.hxx>
30 #include <Standard_Failure.hxx>
31 #include <Standard_OutOfRange.hxx>
32 #include <StdFail_NotDone.hxx>
33 #include <TColgp_Array1OfPnt.hxx>
34 #include <TColStd_Array1OfReal.hxx>
36 //==================================================================================
37 // function : ComputePasInit
38 // purpose : estimate of max step : To avoid abrupt changes during change of isos
39 //==================================================================================
40 void IntWalk_PWalking::ComputePasInit(const Standard_Real theDeltaU1,
41 const Standard_Real theDeltaV1,
42 const Standard_Real theDeltaU2,
43 const Standard_Real theDeltaV2)
45 const Standard_Real aRangePart = 0.01;
46 const Standard_Real Increment = 2.0*pasMax;
47 const Handle(Adaptor3d_HSurface)&
48 Caro1 = myIntersectionOn2S.Function().AuxillarSurface1();
49 const Handle(Adaptor3d_HSurface)&
50 Caro2 = myIntersectionOn2S.Function().AuxillarSurface2();
52 const Standard_Real aDeltaU1=Abs(UM1-Um1);
53 const Standard_Real aDeltaV1=Abs(VM1-Vm1);
54 const Standard_Real aDeltaU2=Abs(UM2-Um2);
55 const Standard_Real aDeltaV2=Abs(VM2-Vm2);
57 //-- limit the reduction of uv box estimate to 0.01 natural box
58 //-- theDeltaU1 : On box of Inter
59 //-- aDeltaU1 : On parametric space
60 if(!Precision::IsInfinite(aDeltaU1))
61 pasuv[0]=Max(Increment*Max(theDeltaU1, aRangePart*aDeltaU1), pasuv[0]);
63 pasuv[0]=Max(Increment*theDeltaU1, pasuv[0]);
65 if(!Precision::IsInfinite(aDeltaV1))
66 pasuv[1]=Max(Increment*Max(theDeltaV1, aRangePart*aDeltaV1), pasuv[1]);
68 pasuv[1]=Max(Increment*theDeltaV1, pasuv[1]);
70 if(!Precision::IsInfinite(aDeltaU2))
71 pasuv[2]=Max(Increment*Max(theDeltaU2, aRangePart*aDeltaU2), pasuv[2]);
73 pasuv[2]=Max(Increment*theDeltaU2, pasuv[2]);
75 if(!Precision::IsInfinite(aDeltaV2))
76 pasuv[3]=Max(Increment*Max(theDeltaV2, aRangePart*aDeltaV2), pasuv[3]);
78 pasuv[3]=Max(Increment*theDeltaV2, pasuv[3]);
80 const Standard_Real ResoU1tol = Adaptor3d_HSurfaceTool::UResolution(Caro1, tolconf);
81 const Standard_Real ResoV1tol = Adaptor3d_HSurfaceTool::VResolution(Caro1, tolconf);
82 const Standard_Real ResoU2tol = Adaptor3d_HSurfaceTool::UResolution(Caro2, tolconf);
83 const Standard_Real ResoV2tol = Adaptor3d_HSurfaceTool::VResolution(Caro2, tolconf);
85 myStepMin[0] = Max(myStepMin[0], 2.0*ResoU1tol);
86 myStepMin[1] = Max(myStepMin[1], 2.0*ResoV1tol);
87 myStepMin[2] = Max(myStepMin[2], 2.0*ResoU2tol);
88 myStepMin[3] = Max(myStepMin[3], 2.0*ResoV2tol);
90 for(Standard_Integer i = 0; i < 4; i++)
92 pasuv[i]=Max(myStepMin[i], pasuv[i]);
96 //=======================================================================
97 //function : IsParallel
98 //purpose : Checks if theLine is parallel of some boundary of given
99 // surface (it is determined by theCheckSurf1 flag).
100 // Parallelism assumes small oscillations (swing is less or
101 // equal than theToler).
102 // Small lines (if first and last parameters in the Surface
103 // are almost equal) are classified as parallel (as same as
104 // any point can be considered as parallel of any line).
105 //=======================================================================
106 static void IsParallel(const Handle(IntSurf_LineOn2S)& theLine,
107 const Standard_Boolean theCheckSurf1,
108 const Standard_Real theToler,
109 Standard_Boolean& theIsUparallel,
110 Standard_Boolean& theIsVparallel)
112 const Standard_Integer aNbPointsMAX = 23;
114 theIsUparallel = theIsVparallel = Standard_True;
116 const Standard_Integer aNbLinePnts = theLine->NbPoints();
117 Standard_Integer aNbPoints = aNbLinePnts;
118 if(aNbPoints > aNbPointsMAX)
120 aNbPoints = aNbPointsMAX;
122 else if(aNbPoints < 3)
124 //Here we cannot estimate parallelism.
125 //Do all same as for small lines
129 Standard_Real aStep = IntToReal(theLine->NbPoints()) / aNbPoints;
130 Standard_Real aNPoint = 1.0;
132 Standard_Real aUmin = RealLast(), aUmax = RealFirst(), aVmin = RealLast(), aVmax = RealFirst();
133 for(Standard_Integer aNum = 1; aNum <= aNbPoints; aNum++, aNPoint += aStep)
135 // Fix possible "out of parameter" case.
136 if (aNPoint > aNbLinePnts)
137 aNPoint = aNbLinePnts;
141 theLine->Value(RealToInt(aNPoint)).ParametersOnS1(u, v);
143 theLine->Value(RealToInt(aNPoint)).ParametersOnS2(u, v);
158 theIsVparallel = ((aUmax - aUmin) < theToler);
159 theIsUparallel = ((aVmax - aVmin) < theToler);
162 //==================================================================================
163 // function : IntWalk_PWalking::IntWalk_PWalking
165 //==================================================================================
166 IntWalk_PWalking::IntWalk_PWalking(const Handle(Adaptor3d_HSurface)& Caro1,
167 const Handle(Adaptor3d_HSurface)& Caro2,
168 const Standard_Real TolTangency,
169 const Standard_Real Epsilon,
170 const Standard_Real Deflection,
171 const Standard_Real Increment )
175 close(Standard_False),
178 myTolTang(TolTangency),
180 myIntersectionOn2S(Caro1,Caro2,TolTangency),
181 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND(0),
182 STATIC_PRECEDENT_INFLEXION(0)
184 Standard_Real KELARG=20.;
186 pasMax=Increment*0.2; //-- June 25 99 after problems with precision
187 Um1 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro1);
188 Vm1 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro1);
189 UM1 = Adaptor3d_HSurfaceTool::LastUParameter(Caro1);
190 VM1 = Adaptor3d_HSurfaceTool::LastVParameter(Caro1);
192 Um2 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro2);
193 Vm2 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro2);
194 UM2 = Adaptor3d_HSurfaceTool::LastUParameter(Caro2);
195 VM2 = Adaptor3d_HSurfaceTool::LastVParameter(Caro2);
197 ResoU1 = Adaptor3d_HSurfaceTool::UResolution(Caro1,Precision::Confusion());
198 ResoV1 = Adaptor3d_HSurfaceTool::VResolution(Caro1,Precision::Confusion());
200 ResoU2 = Adaptor3d_HSurfaceTool::UResolution(Caro2,Precision::Confusion());
201 ResoV2 = Adaptor3d_HSurfaceTool::VResolution(Caro2,Precision::Confusion());
203 Standard_Real NEWRESO;
204 Standard_Real MAXVAL;
205 Standard_Real MAXVAL2;
207 MAXVAL = Abs(Um1); MAXVAL2 = Abs(UM1);
208 if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
209 NEWRESO = ResoU1 * MAXVAL ;
210 if(NEWRESO > ResoU1 &&NEWRESO<10) { ResoU1 = NEWRESO; }
213 MAXVAL = Abs(Um2); MAXVAL2 = Abs(UM2);
214 if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
215 NEWRESO = ResoU2 * MAXVAL ;
216 if(NEWRESO > ResoU2 && NEWRESO<10) { ResoU2 = NEWRESO; }
219 MAXVAL = Abs(Vm1); MAXVAL2 = Abs(VM1);
220 if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
221 NEWRESO = ResoV1 * MAXVAL ;
222 if(NEWRESO > ResoV1 && NEWRESO<10) { ResoV1 = NEWRESO; }
225 MAXVAL = Abs(Vm2); MAXVAL2 = Abs(VM2);
226 if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
227 NEWRESO = ResoV2 * MAXVAL ;
228 if(NEWRESO > ResoV2 && NEWRESO<10) { ResoV2 = NEWRESO; }
230 pasuv[0]=pasMax*Abs(UM1-Um1);
231 pasuv[1]=pasMax*Abs(VM1-Vm1);
232 pasuv[2]=pasMax*Abs(UM2-Um2);
233 pasuv[3]=pasMax*Abs(VM2-Vm2);
235 if(ResoU1>0.0001*pasuv[0]) ResoU1=0.00001*pasuv[0];
236 if(ResoV1>0.0001*pasuv[1]) ResoV1=0.00001*pasuv[1];
237 if(ResoU2>0.0001*pasuv[2]) ResoU2=0.00001*pasuv[2];
238 if(ResoV2>0.0001*pasuv[3]) ResoV2=0.00001*pasuv[3];
241 if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro1)==Standard_False) {
242 //UM1+=KELARG*pasuv[0]; Um1-=KELARG*pasuv[0];
245 Standard_Real t = UM1-Um1;
246 if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro1)) {
247 t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro1)-t);
248 t=(t>KELARG*pasuv[0])? KELARG*pasuv[0] : t;
253 if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro1)==Standard_False) {
254 //VM1+=KELARG*pasuv[1]; Vm1-=KELARG*pasuv[1];
257 Standard_Real t = VM1-Vm1;
258 if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro1)) {
259 t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro1)-t);
260 t=(t>KELARG*pasuv[1])? KELARG*pasuv[1] : t;
265 if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro2)==Standard_False) {
266 //UM2+=KELARG*pasuv[2]; Um2-=KELARG*pasuv[2];
269 Standard_Real t = UM2-Um2;
270 if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro2)) {
271 t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro2)-t);
272 t=(t>KELARG*pasuv[2])? KELARG*pasuv[2] : t;
277 if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro2)==Standard_False) {
278 //VM2+=KELARG*pasuv[3]; Vm2-=KELARG*pasuv[3];
281 Standard_Real t = VM2-Vm2;
282 if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro2)) {
283 t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro2)-t);
284 t=(t>KELARG*pasuv[3])? KELARG*pasuv[3] : t;
289 myStepMin[0] = 100.0*ResoU1;
290 myStepMin[1] = 100.0*ResoV1;
291 myStepMin[2] = 100.0*ResoU2;
292 myStepMin[3] = 100.0*ResoV2;
294 //-- ComputePasInit(pasuv,Um1,UM1,Vm1,VM1,Um2,UM2,Vm2,VM2,Caro1,Caro2);
296 for (Standard_Integer i = 0; i<=3;i++) {
299 pasInit[i] = pasSav[i] = pasuv[i];
304 //==================================================================================
305 // function : IntWalk_PWalking
307 //==================================================================================
308 IntWalk_PWalking::IntWalk_PWalking(const Handle(Adaptor3d_HSurface)& Caro1,
309 const Handle(Adaptor3d_HSurface)& Caro2,
310 const Standard_Real TolTangency,
311 const Standard_Real Epsilon,
312 const Standard_Real Deflection,
313 const Standard_Real Increment,
314 const Standard_Real U1,
315 const Standard_Real V1,
316 const Standard_Real U2,
317 const Standard_Real V2)
321 close(Standard_False),
324 myTolTang(TolTangency),
326 myIntersectionOn2S(Caro1,Caro2,TolTangency),
327 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND(0),
328 STATIC_PRECEDENT_INFLEXION(0)
330 Standard_Real KELARG=20.;
332 pasMax=Increment*0.2; //-- June 25 99 after problems with precision
334 Um1 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro1);
335 Vm1 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro1);
336 UM1 = Adaptor3d_HSurfaceTool::LastUParameter(Caro1);
337 VM1 = Adaptor3d_HSurfaceTool::LastVParameter(Caro1);
339 Um2 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro2);
340 Vm2 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro2);
341 UM2 = Adaptor3d_HSurfaceTool::LastUParameter(Caro2);
342 VM2 = Adaptor3d_HSurfaceTool::LastVParameter(Caro2);
344 ResoU1 = Adaptor3d_HSurfaceTool::UResolution(Caro1,Precision::Confusion());
345 ResoV1 = Adaptor3d_HSurfaceTool::VResolution(Caro1,Precision::Confusion());
347 ResoU2 = Adaptor3d_HSurfaceTool::UResolution(Caro2,Precision::Confusion());
348 ResoV2 = Adaptor3d_HSurfaceTool::VResolution(Caro2,Precision::Confusion());
350 Standard_Real NEWRESO, MAXVAL, MAXVAL2;
354 if(MAXVAL2 > MAXVAL) {
357 NEWRESO = ResoU1 * MAXVAL ;
358 if(NEWRESO > ResoU1) {
364 if(MAXVAL2 > MAXVAL){
367 NEWRESO = ResoU2 * MAXVAL ;
368 if(NEWRESO > ResoU2) {
374 if(MAXVAL2 > MAXVAL) {
377 NEWRESO = ResoV1 * MAXVAL ;
378 if(NEWRESO > ResoV1) {
384 if(MAXVAL2 > MAXVAL){
387 NEWRESO = ResoV2 * MAXVAL ;
388 if(NEWRESO > ResoV2) {
392 pasuv[0]=pasMax*Abs(UM1-Um1);
393 pasuv[1]=pasMax*Abs(VM1-Vm1);
394 pasuv[2]=pasMax*Abs(UM2-Um2);
395 pasuv[3]=pasMax*Abs(VM2-Vm2);
397 if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro1)==Standard_False) {
398 UM1+=KELARG*pasuv[0];
399 Um1-=KELARG*pasuv[0];
402 Standard_Real t = UM1-Um1;
403 if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro1)) {
404 t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro1)-t);
405 t=(t>KELARG*pasuv[0])? KELARG*pasuv[0] : t;
411 if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro1)==Standard_False) {
412 VM1+=KELARG*pasuv[1];
413 Vm1-=KELARG*pasuv[1];
416 Standard_Real t = VM1-Vm1;
417 if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro1)) {
418 t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro1)-t);
419 t=(t>KELARG*pasuv[1])? KELARG*pasuv[1] : t;
424 if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro2)==Standard_False) {
425 UM2+=KELARG*pasuv[2];
426 Um2-=KELARG*pasuv[2];
429 Standard_Real t = UM2-Um2;
430 if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro2)) {
431 t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro2)-t);
432 t=(t>KELARG*pasuv[2])? KELARG*pasuv[2] : t;
438 if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro2)==Standard_False) {
439 VM2+=KELARG*pasuv[3];
440 Vm2-=KELARG*pasuv[3];
443 Standard_Real t = VM2-Vm2;
444 if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro2)) {
445 t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro2)-t);
446 t=(t>KELARG*pasuv[3])? KELARG*pasuv[3] : t;
451 //-- ComputePasInit(pasuv,Um1,UM1,Vm1,VM1,Um2,UM2,Vm2,VM2,Caro1,Caro2);
453 for (Standard_Integer i = 0; i<=3;i++) {
454 pasInit[i] = pasSav[i] = pasuv[i];
457 if(ResoU1>0.0001*pasuv[0]) ResoU1=0.00001*pasuv[0];
458 if(ResoV1>0.0001*pasuv[1]) ResoV1=0.00001*pasuv[1];
459 if(ResoU2>0.0001*pasuv[2]) ResoU2=0.00001*pasuv[2];
460 if(ResoV2>0.0001*pasuv[3]) ResoV2=0.00001*pasuv[3];
462 myStepMin[0] = 100.0*ResoU1;
463 myStepMin[1] = 100.0*ResoV1;
464 myStepMin[2] = 100.0*ResoU2;
465 myStepMin[3] = 100.0*ResoV2;
468 TColStd_Array1OfReal Par(1,4);
476 //==================================================================================
477 // function : PerformFirstPoint
479 //==================================================================================
480 Standard_Boolean IntWalk_PWalking::PerformFirstPoint (const TColStd_Array1OfReal& ParDep,
481 IntSurf_PntOn2S& FirstPoint)
484 close = Standard_False;
487 TColStd_Array1OfReal Param(1,4);
489 for (i=1; i<=4; ++i) {
490 Param(i) = ParDep(i);
492 //-- calculate the first solution point
493 math_FunctionSetRoot Rsnld(myIntersectionOn2S.Function());
495 myIntersectionOn2S.Perform(Param,Rsnld);
496 if (!myIntersectionOn2S.IsDone()) {
497 return Standard_False;
500 if (myIntersectionOn2S.IsEmpty()) {
501 return Standard_False;
504 FirstPoint = myIntersectionOn2S.Point();
505 return Standard_True;
507 //==================================================================================
508 // function : Perform
510 //==================================================================================
511 void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep)
513 Perform(ParDep,Um1,Vm1,Um2,Vm2,UM1,VM1,UM2,VM2);
516 //=======================================================================
517 //function : SQDistPointSurface
518 //purpose : Returns square distance between thePnt and theSurf.
519 // (theU0, theV0) is initial point for extrema
520 //=======================================================================
521 static Standard_Real SQDistPointSurface(const gp_Pnt &thePnt,
522 const Adaptor3d_Surface& theSurf,
523 const Standard_Real theU0,
524 const Standard_Real theV0)
526 Extrema_GenLocateExtPS aExtPS(theSurf);
527 aExtPS.Perform(thePnt, theU0, theV0);
532 return aExtPS.SquareDistance();
535 //==================================================================================
536 // function : IsTangentExtCheck
537 // purpose : Additional check if the surfaces are tangent.
538 // Checks if any point in one surface lie in another surface
539 // (with given tolerance)
540 //==================================================================================
541 static Standard_Boolean IsTangentExtCheck(const Handle(Adaptor3d_HSurface)& theSurf1,
542 const Handle(Adaptor3d_HSurface)& theSurf2,
543 const Standard_Real theU10,
544 const Standard_Real theV10,
545 const Standard_Real theU20,
546 const Standard_Real theV20,
547 const Standard_Real theToler,
548 const Standard_Real theArrStep[])
552 gp_Vec aDu1, aDv1, aDu2, aDv2;
553 theSurf1->D1(theU10, theV10, aPt, aDu1, aDv1);
554 theSurf2->D1(theU20, theV20, aPt, aDu2, aDv2);
556 const gp_Vec aN1(aDu1.Crossed(aDv1)),
557 aN2(aDu2.Crossed(aDv2));
558 const Standard_Real aDP = aN1.Dot(aN2),
559 aSQ1 = aN1.SquareMagnitude(),
560 aSQ2 = aN2.SquareMagnitude();
562 if((aSQ1 < RealSmall()) || (aSQ2 < RealSmall()))
563 return Standard_True; //Tangent
565 if(aDP*aDP < 0.9998*aSQ1*aSQ2)
566 {//cos(ang N1<->N2) < 0.9999
567 return Standard_False; //Not tangent
571 //For two faces (2^2 = 4)
572 const Standard_Real aSQToler = 4.0*theToler*theToler;
573 const Standard_Integer aNbItems = 4;
574 const Standard_Real aParUS1[aNbItems] = { theU10 + theArrStep[0],
575 theU10 - theArrStep[0],
577 const Standard_Real aParVS1[aNbItems] = { theV10, theV10,
578 theV10 + theArrStep[1],
579 theV10 - theArrStep[1]};
580 const Standard_Real aParUS2[aNbItems] = { theU20 + theArrStep[2],
581 theU20 - theArrStep[2],
583 const Standard_Real aParVS2[aNbItems] = { theV20, theV20,
584 theV20 + theArrStep[3],
585 theV20 - theArrStep[3]};
587 for(Standard_Integer i = 0; i < aNbItems; i++)
589 gp_Pnt aP(theSurf1->Value(aParUS1[i], aParVS1[i]));
590 const Standard_Real aSqDist = SQDistPointSurface(aP, theSurf2->Surface(), theU20, theV20);
591 if(aSqDist > aSQToler)
592 return Standard_False;
595 for(Standard_Integer i = 0; i < aNbItems; i++)
597 gp_Pnt aP(theSurf2->Value(aParUS2[i], aParVS2[i]));
598 const Standard_Real aSqDist = SQDistPointSurface(aP, theSurf1->Surface(), theU10, theV10);
599 if(aSqDist > aSQToler)
600 return Standard_False;
603 return Standard_True;
606 //==================================================================================
607 // function : Perform
609 //==================================================================================
610 void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
611 const Standard_Real u1min,
612 const Standard_Real v1min,
613 const Standard_Real u2min,
614 const Standard_Real v2min,
615 const Standard_Real u1max,
616 const Standard_Real v1max,
617 const Standard_Real u2max,
618 const Standard_Real v2max)
620 const Standard_Real aSQDistMax = 1.0e-14;
623 Standard_Integer NbPasOKConseq=0;
624 TColStd_Array1OfReal Param(1,4);
625 IntImp_ConstIsoparametric ChoixIso;
628 done = Standard_False;
631 const Handle(Adaptor3d_HSurface)& Caro1 =myIntersectionOn2S.Function().AuxillarSurface1();
632 const Handle(Adaptor3d_HSurface)& Caro2 =myIntersectionOn2S.Function().AuxillarSurface2();
634 const Standard_Real UFirst1 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro1);
635 const Standard_Real VFirst1 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro1);
636 const Standard_Real ULast1 = Adaptor3d_HSurfaceTool::LastUParameter (Caro1);
637 const Standard_Real VLast1 = Adaptor3d_HSurfaceTool::LastVParameter (Caro1);
639 const Standard_Real UFirst2 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro2);
640 const Standard_Real VFirst2 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro2);
641 const Standard_Real ULast2 = Adaptor3d_HSurfaceTool::LastUParameter (Caro2);
642 const Standard_Real VLast2 = Adaptor3d_HSurfaceTool::LastVParameter (Caro2);
644 ComputePasInit(u1max - u1min,v1max - v1min,u2max - u2min,v2max - v2min);
646 for (Standard_Integer i=0; i<4; ++i)
653 pasInit[i] = pasSav[i] = pasuv[i];
656 line = new IntSurf_LineOn2S ();
658 for (Standard_Integer i=1; i<=4; ++i)
662 //-- reproduce steps uv connected to surfaces Caro1 and Caro2
663 //-- pasuv[] and pasSav[] are modified during the marching
664 for(Standard_Integer i = 0; i < 4; ++i)
666 pasSav[i] = pasuv[i] = pasInit[i];
669 //-- calculate the first solution point
670 math_FunctionSetRoot Rsnld(myIntersectionOn2S.Function());
672 ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld);
673 if (!myIntersectionOn2S.IsDone())
679 if (myIntersectionOn2S.IsEmpty())
684 if(myIntersectionOn2S.IsTangent())
689 Standard_Boolean Arrive, DejaReparti;
690 const Standard_Integer RejectIndexMAX = 250000;
691 Standard_Integer IncKey, RejectIndex;
694 DejaReparti = Standard_False;
698 previousPoint = myIntersectionOn2S.Point();
699 previoustg = Standard_False;
700 previousd = myIntersectionOn2S.Direction();
701 previousd1 = myIntersectionOn2S.DirectionOnS1();
702 previousd2 = myIntersectionOn2S.DirectionOnS2();
705 firstd1 = previousd1;
706 firstd2 = previousd2;
707 tgfirst = tglast = Standard_False;
708 choixIsoSav = ChoixIso;
709 //------------------------------------------------------------
710 //-- Test if the first point of marching corresponds
711 //-- to a point on borders.
712 //-- In this case, DejaReparti is initialized as True
714 pf = previousPoint.Value();
715 Standard_Boolean bTestFirstPoint = Standard_True;
717 previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
719 if(IsTangentExtCheck(Caro1, Caro2, Param(1), Param(2), Param(3), Param(4), myTolTang, pasuv))
722 AddAPoint(line,previousPoint);
724 IntWalk_StatusDeflection Status = IntWalk_OK, aPrevStatus = IntWalk_OK;
725 Standard_Boolean NoTestDeflection = Standard_False;
726 Standard_Real SvParam[4], f;
727 Standard_Integer LevelOfEmptyInmyIntersectionOn2S=0;
728 Standard_Integer LevelOfPointConfondu = 0;
729 Standard_Integer LevelOfIterWithoutAppend = -1;
732 const Standard_Real aTol[4] = { Epsilon(u1max - u1min),
733 Epsilon(v1max - v1min),
734 Epsilon(u2max - u2min),
735 Epsilon(v2max - v2min)};
736 Arrive = Standard_False;
739 aPrevStatus = Status;
741 LevelOfIterWithoutAppend++;
742 if(LevelOfIterWithoutAppend>20)
744 Arrive = Standard_True;
748 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
749 LevelOfIterWithoutAppend = 0;
755 case IntImp_UIsoparametricOnCaro1: f = Abs(previousd1.X()); break;
756 case IntImp_VIsoparametricOnCaro1: f = Abs(previousd1.Y()); break;
757 case IntImp_UIsoparametricOnCaro2: f = Abs(previousd2.X()); break;
758 case IntImp_VIsoparametricOnCaro2: f = Abs(previousd2.Y()); break;
766 previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
769 Standard_Real aIncKey, aEps, dP1, dP2, dP3, dP4;
771 dP1 = sensCheminement * pasuv[0] * previousd1.X() /f;
772 dP2 = sensCheminement * pasuv[1] * previousd1.Y() /f;
773 dP3 = sensCheminement * pasuv[2] * previousd2.X() /f;
774 dP4 = sensCheminement * pasuv[3] * previousd2.Y() /f;
776 aIncKey=5.*(Standard_Real)IncKey;
778 if(ChoixIso == IntImp_UIsoparametricOnCaro1 && Abs(dP1) < aEps)
783 if(ChoixIso == IntImp_VIsoparametricOnCaro1 && Abs(dP2) < aEps)
788 if(ChoixIso == IntImp_UIsoparametricOnCaro2 && Abs(dP3) < aEps)
793 if(ChoixIso == IntImp_VIsoparametricOnCaro2 && Abs(dP4) < aEps)
803 //==========================
809 Standard_Integer aTryNumber = 0;
810 Standard_Real isBadPoint = Standard_False;
811 IntImp_ConstIsoparametric aBestIso = ChoixIso;
814 isBadPoint = Standard_False;
816 ChoixIso= myIntersectionOn2S.Perform(Param, Rsnld, aBestIso);
818 if (myIntersectionOn2S.IsDone() && !myIntersectionOn2S.IsEmpty())
820 //If we go along any surface boundary then it is possible
821 //to find "outboundaried" point.
822 //Nevertheless, if this deflection is quite small, we will be
823 //able to adjust this point to the boundary.
825 Standard_Real aNewPnt[4], anAbsParamDist[4];
826 myIntersectionOn2S.Point().Parameters(aNewPnt[0], aNewPnt[1], aNewPnt[2], aNewPnt[3]);
827 const Standard_Real aParMin[4] = {u1min, v1min, u2min, v2min};
828 const Standard_Real aParMax[4] = {u1max, v1max, u2max, v2max};
830 for(Standard_Integer i = 0; i < 4; i++)
832 if(Abs(aNewPnt[i] - aParMin[i]) < aTol[i])
833 aNewPnt[i] = aParMin[i];
834 else if(Abs(aNewPnt[i] - aParMax[i]) < aTol[i])
835 aNewPnt[i] = aParMax[i];
838 if (aNewPnt[0] < u1min || aNewPnt[0] > u1max ||
839 aNewPnt[1] < v1min || aNewPnt[1] > v1max ||
840 aNewPnt[2] < u2min || aNewPnt[2] > u2max ||
841 aNewPnt[3] < v2min || aNewPnt[3] > v2max)
843 break; // Out of borders, handle this later.
846 myIntersectionOn2S.ChangePoint().SetValue(aNewPnt[0],
851 anAbsParamDist[0] = Abs(Param(1) - dP1 - aNewPnt[0]);
852 anAbsParamDist[1] = Abs(Param(2) - dP2 - aNewPnt[1]);
853 anAbsParamDist[2] = Abs(Param(3) - dP3 - aNewPnt[2]);
854 anAbsParamDist[3] = Abs(Param(4) - dP4 - aNewPnt[3]);
855 if (anAbsParamDist[0] < ResoU1 &&
856 anAbsParamDist[1] < ResoV1 &&
857 anAbsParamDist[2] < ResoU2 &&
858 anAbsParamDist[3] < ResoV2 &&
859 Status != IntWalk_PasTropGrand)
861 isBadPoint = Standard_True;
862 aBestIso = IntImp_ConstIsoparametric((aBestIso + 1) % 4);
865 } while (isBadPoint && ++aTryNumber <= 4);
867 if (!myIntersectionOn2S.IsDone())
869 //end of line, division
870 Arrive = Standard_False;
875 RepartirOuDiviser(DejaReparti, ChoixIso, Arrive);
879 //== Calculation of exact point from Param(.) is possible
880 if (myIntersectionOn2S.IsEmpty())
882 Standard_Real u1,v1,u2,v2;
883 previousPoint.Parameters(u1,v1,u2,v2);
885 Arrive = Standard_False;
886 if(u1<UFirst1 || u1>ULast1)
888 Arrive=Standard_True;
891 if(u2<UFirst2 || u2>ULast2)
893 Arrive=Standard_True;
896 if(v1<VFirst1 || v1>VLast1)
898 Arrive=Standard_True;
901 if(v2<VFirst2 || v2>VLast2)
903 Arrive=Standard_True;
906 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
907 LevelOfEmptyInmyIntersectionOn2S++;
909 if(LevelOfEmptyInmyIntersectionOn2S>10)
919 //============================================================
920 //== A point has been found : T E S T D E F L E C T I O N
921 //============================================================
924 NoTestDeflection = Standard_False;
928 if(--LevelOfEmptyInmyIntersectionOn2S<=0)
930 LevelOfEmptyInmyIntersectionOn2S=0;
931 if(LevelOfIterWithoutAppend < 10)
933 Status = TestDeflection(ChoixIso);
945 //============================================================
946 //== T r a i t e m e n t s u r S t a t u s ==
947 //============================================================
948 if(LevelOfPointConfondu > 5)
950 Status = IntWalk_ArretSurPoint;
951 LevelOfPointConfondu = 0;
954 if(Status==IntWalk_OK)
957 if(NbPasOKConseq >= 5)
960 Standard_Boolean pastroppetit;
965 pastroppetit=Standard_True;
967 if(pasuv[0]<pasInit[0])
969 t = (pasInit[0]-pasuv[0])*0.25;
976 pastroppetit=Standard_False;
979 if(pasuv[1]<pasInit[1])
981 t = (pasInit[1]-pasuv[1])*0.25;
982 if(t>0.1*pasInit[1]) {
987 pastroppetit=Standard_False;
990 if(pasuv[2]<pasInit[2])
992 t = (pasInit[2]-pasuv[2])*0.25;
999 pastroppetit=Standard_False;
1002 if(pasuv[3]<pasInit[3])
1004 t = (pasInit[3]-pasuv[3])*0.25;
1005 if(t>0.1*pasInit[3]) {
1009 pastroppetit=Standard_False;
1023 pastroppetit=Standard_False;
1027 while(pastroppetit);
1029 }//Status==IntWalk_OK
1036 case IntWalk_ArretSurPointPrecedent:
1038 Arrive = Standard_False;
1039 RepartirOuDiviser(DejaReparti, ChoixIso, Arrive);
1042 case IntWalk_PasTropGrand:
1044 Param(1)=SvParam[0];
1045 Param(2)=SvParam[1];
1046 Param(3)=SvParam[2];
1047 Param(4)=SvParam[3];
1049 if(LevelOfIterWithoutAppend > 5)
1051 for (Standard_Integer i = 0; i < 4; i++)
1053 if (pasSav[i] > pasInit[i])
1056 const Standard_Real aDelta = (pasInit[i]-pasSav[i])*0.25;
1058 if(aDelta > Epsilon(pasInit[i]))
1060 pasInit[i] -= aDelta;
1061 LevelOfIterWithoutAppend=0;
1068 case IntWalk_PointConfondu:
1070 LevelOfPointConfondu++;
1072 if(LevelOfPointConfondu>5)
1074 Standard_Boolean pastroppetit;
1078 pastroppetit=Standard_True;
1080 for(Standard_Integer i = 0; i < 4; i++)
1082 if(pasuv[i]<pasInit[i])
1084 pasuv[i]+=(pasInit[i]-pasuv[i])*0.25;
1085 pastroppetit=Standard_False;
1101 pastroppetit=Standard_False;
1105 while(pastroppetit);
1110 case IntWalk_StepTooSmall:
1112 Standard_Boolean hasStepBeenIncreased = Standard_False;
1114 for(Standard_Integer i = 0; i < 4; i++)
1116 const Standard_Real aNewStep = Min(1.5*pasuv[i], pasInit[i]);
1117 if(aNewStep > pasuv[i])
1119 pasuv[i] = aNewStep;
1120 hasStepBeenIncreased = Standard_True;
1124 if(hasStepBeenIncreased)
1126 Param(1)=SvParam[0];
1127 Param(2)=SvParam[1];
1128 Param(3)=SvParam[2];
1129 Param(4)=SvParam[3];
1131 // In order to avoid cyclic changes
1132 // (PasTropGrand --> Decrease step -->
1133 // StepTooSmall --> Increase step --> PasTropGrand...)
1134 // nullify LevelOfIterWithoutAppend only if the condition
1136 if (aPrevStatus != IntWalk_PasTropGrand)
1137 LevelOfIterWithoutAppend = 0;
1143 case IntWalk_ArretSurPoint://006
1145 //=======================================================
1146 //== Stop Test t : Frame on Param(.) ==
1147 //=======================================================
1149 Arrive = TestArret(DejaReparti,Param,ChoixIso);
1150 // JMB 30th December 1999.
1151 // Some statement below should not be put in comment because they are useful.
1152 // See grid CTO 909 A1 which infinitely loops
1153 if(Arrive==Standard_False && Status==IntWalk_ArretSurPoint)
1155 Arrive=Standard_True;
1157 cout << "IntWalk_PWalking_1.gxx: Problems with intersection"<<endl;
1163 NbPasOKConseq = -10;
1168 //=====================================================
1169 //== Param(.) is in the limits ==
1170 //== and does not end a closed line ==
1171 //=====================================================
1172 //== Check on the current point of myInters
1173 Standard_Boolean pointisvalid = Standard_False;
1175 Standard_Real u1,v1,u2,v2;
1176 myIntersectionOn2S.Point().Parameters(u1,v1,u2,v2);
1179 if(u1 <= UM1 && u2 <= UM2 && v1 <= VM1 &&
1180 v2 <= VM2 && u1 >= Um1 && u2 >= Um2 &&
1181 v1 >= Vm1 && v2 >= Vm2)
1183 pointisvalid=Standard_True;
1190 previousPoint = myIntersectionOn2S.Point();
1191 previoustg = myIntersectionOn2S.IsTangent();
1195 previousd = myIntersectionOn2S.Direction();
1196 previousd1 = myIntersectionOn2S.DirectionOnS1();
1197 previousd2 = myIntersectionOn2S.DirectionOnS2();
1199 //=====================================================
1200 //== Check on the previous Point
1202 Standard_Real u1,v1,u2,v2;
1203 previousPoint.Parameters(u1,v1,u2,v2);
1204 if( u1 <= UM1 && u2 <= UM2 && v1 <= VM1 &&
1205 v2 <= VM2 && u1 >= Um1 && u2 >= Um2 &&
1206 v1 >= Vm1 && v2 >= Vm2)
1208 pl = previousPoint.Value();
1211 if(pf.SquareDistance(pl) < aSQDistMax)
1221 bTestFirstPoint = Standard_False;
1225 AddAPoint(line,previousPoint);
1228 if(RejectIndex >= RejectIndexMAX)
1230 Arrive = Standard_True;
1235 LevelOfIterWithoutAppend = 0;
1239 //====================================================
1241 if(Status == IntWalk_ArretSurPoint)
1243 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1247 if (line->NbPoints() == 2)
1249 pasSav[0] = pasuv[0];
1250 pasSav[1] = pasuv[1];
1251 pasSav[2] = pasuv[2];
1252 pasSav[3] = pasuv[3];
1260 //================= la ligne est fermee ===============
1261 AddAPoint(line,line->Value(1)); //ligne fermee
1262 LevelOfIterWithoutAppend=0;
1266 //====================================================
1267 //== Param was not in the limits (was reframed)
1268 //====================================================
1269 Standard_Boolean bPrevNotTangent = !previoustg || !myIntersectionOn2S.IsTangent();
1271 IntImp_ConstIsoparametric SauvChoixIso = ChoixIso;
1272 ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld,ChoixIso);
1274 if(!myIntersectionOn2S.IsEmpty()) //002
1276 // mutially outpasses in the square or intersection in corner
1278 if(TestArret(Standard_True,Param,ChoixIso))
1280 NbPasOKConseq = -10;
1281 ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld,ChoixIso);
1283 if(!myIntersectionOn2S.IsEmpty())
1285 previousPoint = myIntersectionOn2S.Point();
1286 previoustg = myIntersectionOn2S.IsTangent();
1290 previousd = myIntersectionOn2S.Direction();
1291 previousd1 = myIntersectionOn2S.DirectionOnS1();
1292 previousd2 = myIntersectionOn2S.DirectionOnS2();
1295 pl = previousPoint.Value();
1299 if(pf.SquareDistance(pl) < aSQDistMax)
1309 bTestFirstPoint = Standard_False;
1313 AddAPoint(line,previousPoint);
1316 if(RejectIndex >= RejectIndexMAX)
1318 Arrive = Standard_True;
1323 LevelOfIterWithoutAppend=0;
1324 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1328 //fail framing divides the step
1329 Arrive = Standard_False;
1330 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1331 NoTestDeflection = Standard_True;
1332 ChoixIso = SauvChoixIso;
1337 // save the last point
1338 // to revert to it if the current point is out of bounds
1340 IntSurf_PntOn2S previousPointSave = previousPoint;
1341 Standard_Boolean previoustgSave = previoustg;
1342 gp_Dir previousdSave = previousd;
1343 gp_Dir2d previousd1Save = previousd1;
1344 gp_Dir2d previousd2Save = previousd2;
1346 previousPoint = myIntersectionOn2S.Point();
1347 previoustg = myIntersectionOn2S.IsTangent();
1348 Arrive = Standard_False;
1352 previousd = myIntersectionOn2S.Direction();
1353 previousd1 = myIntersectionOn2S.DirectionOnS1();
1354 previousd2 = myIntersectionOn2S.DirectionOnS2();
1357 //========================================
1358 //== Check on PreviousPoint @@
1361 Standard_Real u1,v1,u2,v2;
1362 previousPoint.Parameters(u1,v1,u2,v2);
1364 //To save initial 2d points
1365 gp_Pnt2d ParamPntOnS1(Param(1), Param(2));
1366 gp_Pnt2d ParamPntOnS2(Param(3), Param(4));
1368 ///////////////////////////
1376 Standard_Boolean bFlag1, bFlag2;
1377 Standard_Real aTol2D=1.e-11;
1379 bFlag1=u1 >= Um1-aTol2D && v1 >= Vm1-aTol2D && u1 <= UM1+aTol2D && v1 <= VM1+aTol2D;
1380 bFlag2=u2 >= Um2-aTol2D && v2 >= Vm2-aTol2D && u2 <= UM2+aTol2D && v2 <= VM2+aTol2D;
1381 if (bFlag1 && bFlag2)
1383 if (line->NbPoints() > 1)
1385 IntSurf_PntOn2S prevprevPoint = line->Value(line->NbPoints()-1);
1386 Standard_Real ppU1, ppV1, ppU2, ppV2;
1387 prevprevPoint.Parameters(ppU1, ppV1, ppU2, ppV2);
1388 Standard_Real pU1, pV1, pU2, pV2;
1389 previousPointSave.Parameters(pU1, pV1, pU2, pV2);
1390 gp_Vec2d V1onS1(gp_Pnt2d(ppU1, ppV1), gp_Pnt2d(pU1, pV1));
1391 gp_Vec2d V2onS1(gp_Pnt2d(pU1, pV1), gp_Pnt2d(u1, v1));
1392 gp_Vec2d V1onS2(gp_Pnt2d(ppU2, ppV2), gp_Pnt2d(pU2, pV2));
1393 gp_Vec2d V2onS2(gp_Pnt2d(pU2, pV2), gp_Pnt2d(u2, v2));
1395 const Standard_Real aDot1 = V1onS1 * V2onS1;
1396 const Standard_Real aDot2 = V1onS2 * V2onS2;
1398 if ((aDot1 < 0.0) || (aDot2 < 0.0))
1400 Arrive = Standard_True;
1405 if(u1 <= UM1 && u2 <= UM2 && v1 <= VM1 &&
1406 v2 <= VM2 && u1 >= Um1 && u2 >= Um2 &&
1407 v1 >= Vm1 && v2 >= Vm2) {
1410 pl = previousPoint.Value();
1414 if(pf.SquareDistance(pl) < aSQDistMax)
1425 bTestFirstPoint = Standard_False;
1429 //To avoid walking around the same point
1430 //in the tangent zone near a border
1434 //There are three consecutive points:
1435 //previousPointSave -> ParamPnt -> curPnt.
1437 Standard_Real prevU1, prevV1, prevU2, prevV2;
1438 previousPointSave.Parameters(prevU1, prevV1, prevU2, prevV2);
1439 gp_Pnt2d prevPntOnS1(prevU1, prevV1), prevPntOnS2(prevU2, prevV2);
1440 gp_Pnt2d curPntOnS1(u1, v1), curPntOnS2(u2, v2);
1441 gp_Vec2d PrevToParamOnS1(prevPntOnS1, ParamPntOnS1);
1442 gp_Vec2d PrevToCurOnS1(prevPntOnS1, curPntOnS1);
1443 gp_Vec2d PrevToParamOnS2(prevPntOnS2, ParamPntOnS2);
1444 gp_Vec2d PrevToCurOnS2(prevPntOnS2, curPntOnS2);
1445 Standard_Real MaxAngle = 3*M_PI/4;
1446 Standard_Real anAngleS1 = 0.0, anAngleS2 = 0.0;
1447 const Standard_Real aSQMParS1 = PrevToParamOnS1.SquareMagnitude();
1448 const Standard_Real aSQMParS2 = PrevToParamOnS2.SquareMagnitude();
1449 const Standard_Real aSQMCurS1 = PrevToCurOnS1.SquareMagnitude();
1450 const Standard_Real aSQMCurS2 = PrevToCurOnS2.SquareMagnitude();
1452 if(aSQMCurS1 < gp::Resolution())
1454 //We came back to the one of previos point.
1455 //Therefore, we must break;
1459 else if(aSQMParS1 < gp::Resolution())
1461 //We are walking along tangent zone.
1462 //It should be continued.
1467 anAngleS1 = Abs(PrevToParamOnS1.Angle(PrevToCurOnS1));
1470 if(aSQMCurS2 < gp::Resolution())
1472 //We came back to the one of previos point.
1473 //Therefore, we must break;
1477 else if(aSQMParS2 < gp::Resolution())
1479 //We are walking along tangent zone.
1480 //It should be continued;
1485 anAngleS2 = Abs(PrevToParamOnS2.Angle(PrevToCurOnS2));
1488 if ((anAngleS1 > MaxAngle) && (anAngleS2 > MaxAngle))
1490 Arrive = Standard_True;
1495 //Check singularity.
1496 //I.e. check if we are walking along direction, which does not
1497 //result in comming to any point (i.e. derivative
1498 //3D-intersection curve along this direction is equal to 0).
1499 //A sphere with direction {dU=1, dV=0} from point
1500 //(U=0, V=M_PI/2) can be considered as example for
1501 //this case (we cannot find another 3D-point if we go thus).
1503 //Direction chosen along 1st and 2nd surface correspondingly
1504 const gp_Vec2d aDirS1(prevPntOnS1, curPntOnS1),
1505 aDirS2(prevPntOnS2, curPntOnS2);
1508 gp_Vec aDuS1, aDvS1, aDuS2, aDvS2;
1510 myIntersectionOn2S.Function().AuxillarSurface1()->
1511 D1(curPntOnS1.X(), curPntOnS1.Y(), aPtemp, aDuS1, aDvS1);
1512 myIntersectionOn2S.Function().AuxillarSurface2()->
1513 D1(curPntOnS2.X(), curPntOnS2.Y(), aPtemp, aDuS2, aDvS2);
1515 //Derivative WLine along (it is vector-function indeed)
1517 //(https://en.wikipedia.org/wiki/Directional_derivative#Variation_using_only_direction_of_vector).
1518 //F1 - on the 1st surface, F2 - on the 2nd surface.
1519 //x, y, z - coordinates of derivative vector.
1520 const Standard_Real aF1x = aDuS1.X()*aDirS1.X() +
1521 aDvS1.X()*aDirS1.Y();
1522 const Standard_Real aF1y = aDuS1.Y()*aDirS1.X() +
1523 aDvS1.Y()*aDirS1.Y();
1524 const Standard_Real aF1z = aDuS1.Z()*aDirS1.X() +
1525 aDvS1.Z()*aDirS1.Y();
1526 const Standard_Real aF2x = aDuS2.X()*aDirS2.X() +
1527 aDvS2.X()*aDirS2.Y();
1528 const Standard_Real aF2y = aDuS2.Y()*aDirS2.X() +
1529 aDvS2.Y()*aDirS2.Y();
1530 const Standard_Real aF2z = aDuS2.Z()*aDirS2.X() +
1531 aDvS2.Z()*aDirS2.Y();
1533 const Standard_Real aF1 = aF1x*aF1x + aF1y*aF1y + aF1z*aF1z;
1534 const Standard_Real aF2 = aF2x*aF2x + aF2y*aF2y + aF2z*aF2z;
1536 if((aF1 < gp::Resolution()) && (aF2 < gp::Resolution()))
1538 //All derivative are equal to 0. Therefore, there is
1539 //no point in going along direction chosen.
1540 Arrive = Standard_True;
1544 }//if (previoustg) cond.
1546 ////////////////////////////////////////
1547 AddAPoint(line,previousPoint);
1550 if(RejectIndex >= RejectIndexMAX)
1552 Arrive = Standard_True;
1558 LevelOfIterWithoutAppend=0;
1559 Arrive = Standard_True;
1563 // revert to the last correctly calculated point
1564 previousPoint = previousPointSave;
1565 previoustg = previoustgSave;
1566 previousd = previousdSave;
1567 previousd1 = previousd1Save;
1568 previousd2 = previousd2Save;
1573 Standard_Boolean wasExtended = Standard_False;
1575 if(Arrive && myIntersectionOn2S.IsTangent() && bPrevNotTangent)
1577 if(ExtendLineInCommonZone(SauvChoixIso, DejaReparti))
1579 wasExtended = Standard_True;
1580 Arrive = Standard_False;
1581 ChoixIso = SauvChoixIso;
1585 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1588 myIntersectionOn2S.IsDone() && !myIntersectionOn2S.IsEmpty() &&
1589 myIntersectionOn2S.IsTangent() && bPrevNotTangent &&
1592 if(ExtendLineInCommonZone(SauvChoixIso, DejaReparti))
1594 wasExtended = Standard_True;
1595 Arrive = Standard_False;
1596 ChoixIso = SauvChoixIso;
1599 }//else !TestArret() $
1600 }//$$ end successful framing on border (!myIntersectionOn2S.IsEmpty())
1603 //echec framing on border; division of step
1604 Arrive = Standard_False;
1605 NoTestDeflection = Standard_True;
1606 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1608 }//$$$ end framing on border (!close)
1609 }//004 fin TestArret return Arrive = True
1610 } // 006case IntWalk_ArretSurPoint: end Processing Status = OK or ArretSurPoint
1611 } //007 switch(Status)
1612 } //008 end processing point (TEST DEFLECTION)
1613 } //009 end processing line (else if myIntersectionOn2S.IsDone())
1614 } //010 end if first departure point allows marching while (!Arrive)
1616 done = Standard_True;
1618 // ===========================================================================================================
1619 // function: ExtendLineInCommonZone
1620 // purpose: Extends already computed line inside tangent zone in the direction given by theChoixIso.
1621 // Returns Standard_True if the line was extended through tangent zone and the last computed point
1622 // is outside the tangent zone (but it is not put into the line). Otherwise returns Standard_False.
1623 // ===========================================================================================================
1624 Standard_Boolean IntWalk_PWalking::ExtendLineInCommonZone(const IntImp_ConstIsoparametric theChoixIso,
1625 const Standard_Boolean theDirectionFlag)
1627 Standard_Boolean bOutOfTangentZone = Standard_False;
1628 Standard_Boolean bStop = !myIntersectionOn2S.IsTangent();
1629 Standard_Integer dIncKey = 1;
1630 TColStd_Array1OfReal Param(1,4);
1631 IntWalk_StatusDeflection Status = IntWalk_OK;
1632 Standard_Integer nbIterWithoutAppend = 0;
1633 Standard_Integer nbEqualPoints = 0;
1634 Standard_Integer parit = 0;
1635 Standard_Integer uvit = 0;
1636 IntSurf_SequenceOfPntOn2S aSeqOfNewPoint;
1639 nbIterWithoutAppend++;
1641 if((nbIterWithoutAppend > 20) || (nbEqualPoints > 20)) {
1643 cout<<"Infinite loop detected. Stop iterations (IntWalk_PWalking_1.gxx)" << endl;
1645 bStop = Standard_True;
1648 Standard_Real f = 0.;
1650 switch (theChoixIso)
1652 case IntImp_UIsoparametricOnCaro1: f = Abs(previousd1.X()); break;
1653 case IntImp_VIsoparametricOnCaro1: f = Abs(previousd1.Y()); break;
1654 case IntImp_UIsoparametricOnCaro2: f = Abs(previousd2.X()); break;
1655 case IntImp_VIsoparametricOnCaro2: f = Abs(previousd2.Y()); break;
1660 previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
1662 Standard_Real dP1 = sensCheminement * pasuv[0] * previousd1.X() /f;
1663 Standard_Real dP2 = sensCheminement * pasuv[1] * previousd1.Y() /f;
1664 Standard_Real dP3 = sensCheminement * pasuv[2] * previousd2.X() /f;
1665 Standard_Real dP4 = sensCheminement * pasuv[3] * previousd2.Y() /f;
1667 if(theChoixIso == IntImp_UIsoparametricOnCaro1 && Abs(dP1) < 1.e-7) dP1 *= (5. * (Standard_Real)dIncKey);
1668 if(theChoixIso == IntImp_VIsoparametricOnCaro1 && Abs(dP2) < 1.e-7) dP2 *= (5. * (Standard_Real)dIncKey);
1669 if(theChoixIso == IntImp_UIsoparametricOnCaro2 && Abs(dP3) < 1.e-7) dP3 *= (5. * (Standard_Real)dIncKey);
1670 if(theChoixIso == IntImp_VIsoparametricOnCaro2 && Abs(dP4) < 1.e-7) dP4 *= (5. * (Standard_Real)dIncKey);
1676 Standard_Real SvParam[4];
1677 IntImp_ConstIsoparametric ChoixIso = theChoixIso;
1679 for(parit = 0; parit < 4; parit++) {
1680 SvParam[parit] = Param(parit+1);
1682 math_FunctionSetRoot Rsnld(myIntersectionOn2S.Function());
1683 ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld, theChoixIso);
1685 if (!myIntersectionOn2S.IsDone()) {
1686 return bOutOfTangentZone;
1689 if (myIntersectionOn2S.IsEmpty()) {
1690 return bOutOfTangentZone;
1693 Status = TestDeflection(ChoixIso);
1695 if(Status == IntWalk_OK) {
1697 for(uvit = 0; uvit < 4; uvit++) {
1698 if(pasuv[uvit] < pasInit[uvit]) {
1699 pasuv[uvit] = pasInit[uvit];
1705 case IntWalk_ArretSurPointPrecedent:
1707 bStop = Standard_True;
1708 bOutOfTangentZone = !myIntersectionOn2S.IsTangent();
1711 case IntWalk_PasTropGrand:
1713 for(parit = 0; parit < 4; parit++) {
1714 Param(parit+1) = SvParam[parit];
1716 Standard_Boolean bDecrease = Standard_False;
1718 for(uvit = 0; uvit < 4; uvit++) {
1719 if(pasSav[uvit] < pasInit[uvit]) {
1720 pasInit[uvit] -= (pasInit[uvit] - pasSav[uvit]) * 0.1;
1721 bDecrease = Standard_True;
1725 if(bDecrease) nbIterWithoutAppend--;
1728 case IntWalk_PointConfondu:
1730 for(uvit = 0; uvit < 4; uvit++) {
1731 if(pasuv[uvit] < pasInit[uvit]) {
1732 pasuv[uvit] += (pasInit[uvit] - pasuv[uvit]) * 0.1;
1738 case IntWalk_ArretSurPoint:
1741 bStop = TestArret(theDirectionFlag, Param, ChoixIso);
1746 Standard_Real u11,v11,u12,v12;
1747 myIntersectionOn2S.Point().Parameters(u11,v11,u12,v12);
1748 Standard_Real u21,v21,u22,v22;
1749 previousPoint.Parameters(u21,v21,u22,v22);
1751 if(((fabs(u11-u21) < ResoU1) && (fabs(v11-v21) < ResoV1)) ||
1752 ((fabs(u12-u22) < ResoU2) && (fabs(v12-v22) < ResoV2))) {
1761 bStop = bStop || !myIntersectionOn2S.IsTangent();
1762 bOutOfTangentZone = !myIntersectionOn2S.IsTangent();
1765 Standard_Boolean pointisvalid = Standard_False;
1766 Standard_Real u1,v1,u2,v2;
1767 myIntersectionOn2S.Point().Parameters(u1,v1,u2,v2);
1769 if(u1 <= UM1 && u2 <= UM2 && v1 <= VM1 &&
1770 v2 <= VM2 && u1 >= Um1 && u2 >= Um2 &&
1771 v1 >= Vm1 && v2 >= Vm2)
1772 pointisvalid = Standard_True;
1775 previousPoint = myIntersectionOn2S.Point();
1776 previoustg = myIntersectionOn2S.IsTangent();
1779 previousd = myIntersectionOn2S.Direction();
1780 previousd1 = myIntersectionOn2S.DirectionOnS1();
1781 previousd2 = myIntersectionOn2S.DirectionOnS2();
1783 Standard_Boolean bAddPoint = Standard_True;
1785 if(line->NbPoints() >= 1) {
1786 gp_Pnt pf = line->Value(1).Value();
1787 gp_Pnt pl = previousPoint.Value();
1789 if(pf.Distance(pl) < Precision::Confusion()) {
1791 if(dIncKey == 5000) return bOutOfTangentZone;
1792 else bAddPoint = Standard_False;
1797 aSeqOfNewPoint.Append(previousPoint);
1798 nbIterWithoutAppend = 0;
1802 if (line->NbPoints() == 2) {
1803 for(uvit = 0; uvit < 4; uvit++) {
1804 pasSav[uvit] = pasuv[uvit];
1808 if ( !pointisvalid ) {
1809 // decrease step if out of bounds
1810 // otherwise the same calculations will be
1811 // repeated several times
1812 if ( ( u1 > UM1 ) || ( u1 < Um1 ) )
1815 if ( ( v1 > VM1 ) || ( v1 < Vm1 ) )
1818 if ( ( u2 > UM2 ) || ( u2 < Um2 ) )
1821 if ( ( v2 > VM2 ) || ( v2 < Vm2 ) )
1826 if(close && (line->NbPoints() >= 1)) {
1828 if(!bOutOfTangentZone) {
1829 aSeqOfNewPoint.Append(line->Value(1)); // line end
1831 nbIterWithoutAppend = 0;
1834 ChoixIso = myIntersectionOn2S.Perform(Param, Rsnld, theChoixIso);
1836 if(myIntersectionOn2S.IsEmpty()) {
1837 bStop = Standard_True;// !myIntersectionOn2S.IsTangent();
1838 bOutOfTangentZone = Standard_False; // !myIntersectionOn2S.IsTangent();
1841 Standard_Boolean bAddPoint = Standard_True;
1842 Standard_Boolean pointisvalid = Standard_False;
1844 previousPoint = myIntersectionOn2S.Point();
1845 Standard_Real u1,v1,u2,v2;
1846 previousPoint.Parameters(u1,v1,u2,v2);
1848 if(u1 <= UM1 && u2 <= UM2 && v1 <= VM1 &&
1849 v2 <= VM2 && u1 >= Um1 && u2 >= Um2 &&
1850 v1 >= Vm1 && v2 >= Vm2)
1851 pointisvalid = Standard_True;
1855 if(line->NbPoints() >= 1) {
1856 gp_Pnt pf = line->Value(1).Value();
1857 gp_Pnt pl = previousPoint.Value();
1859 if(pf.Distance(pl) < Precision::Confusion()) {
1861 if(dIncKey == 5000) return bOutOfTangentZone;
1862 else bAddPoint = Standard_False;
1866 if(bAddPoint && !bOutOfTangentZone) {
1867 aSeqOfNewPoint.Append(previousPoint);
1868 nbIterWithoutAppend = 0;
1883 Standard_Boolean bExtendLine = Standard_False;
1884 Standard_Real u1 = 0., v1 = 0., u2 = 0., v2 = 0.;
1886 Standard_Integer pit = 0;
1888 for(pit = 0; !bExtendLine && (pit < 2); pit++) {
1890 previousPoint.Parameters(u1,v1,u2,v2);
1892 if(aSeqOfNewPoint.Length() > 0)
1893 aSeqOfNewPoint.Value(aSeqOfNewPoint.Length()).Parameters(u1,v1,u2,v2);
1898 if(((u1 - Um1) < ResoU1) ||
1899 ((UM1 - u1) < ResoU1) ||
1900 ((u2 - Um2) < ResoU2) ||
1901 ((UM2 - u2) < ResoU2) ||
1902 ((v1 - Vm1) < ResoV1) ||
1903 ((VM1 - v1) < ResoV1) ||
1904 ((v2 - Vm2) < ResoV2) ||
1905 ((VM2 - v2) < ResoV2))
1906 bExtendLine = Standard_True;
1910 // if(Status == IntWalk_OK || Status == IntWalk_ArretSurPoint) {
1911 if(Status == IntWalk_OK) {
1912 bExtendLine = Standard_True;
1914 if(aSeqOfNewPoint.Length() > 1) {
1915 TColStd_Array1OfReal FirstParams(0, 3), LastParams(0, 3), Resolutions(0, 3);
1916 Resolutions(0) = ResoU1; Resolutions(1) = ResoV1; Resolutions(2) = ResoU2; Resolutions(3) = ResoV2;
1918 aSeqOfNewPoint(1).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
1919 FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
1920 aSeqOfNewPoint(aSeqOfNewPoint.Length()).Parameters(LastParams.ChangeValue(0),
1921 LastParams.ChangeValue(1),
1922 LastParams.ChangeValue(2),
1923 LastParams.ChangeValue(3));
1924 Standard_Integer indexofiso = 0;
1926 if(theChoixIso == IntImp_UIsoparametricOnCaro1) indexofiso = 0;
1927 if(theChoixIso == IntImp_VIsoparametricOnCaro1) indexofiso = 1;
1928 if(theChoixIso == IntImp_UIsoparametricOnCaro2) indexofiso = 2;
1929 if(theChoixIso == IntImp_VIsoparametricOnCaro2) indexofiso = 3;
1931 Standard_Integer afirstindex = (indexofiso < 2) ? 0 : 2;
1932 gp_Vec2d aTangentZoneDir(gp_Pnt2d(FirstParams.Value(afirstindex), FirstParams.Value(afirstindex + 1)),
1933 gp_Pnt2d(LastParams.Value(afirstindex), LastParams.Value(afirstindex + 1)));
1935 gp_Dir2d anIsoDir(0, 1);
1937 if((indexofiso == 1) || (indexofiso == 3))
1938 anIsoDir = gp_Dir2d(1, 0);
1940 if(aTangentZoneDir.SquareMagnitude() > gp::Resolution()) {
1941 Standard_Real piquota = M_PI*0.25;
1943 if(fabs(aTangentZoneDir.Angle(anIsoDir)) > piquota) {
1944 Standard_Integer ii = 1, nextii = 2;
1946 Standard_Real asqresol = gp::Resolution();
1947 asqresol *= asqresol;
1950 aSeqOfNewPoint(ii).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
1951 FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
1952 aSeqOfNewPoint(ii + 1).Parameters(LastParams.ChangeValue(0), LastParams.ChangeValue(1),
1953 LastParams.ChangeValue(2), LastParams.ChangeValue(3));
1954 d1 = gp_Vec2d(gp_Pnt2d(FirstParams.Value(afirstindex),
1955 FirstParams.Value(afirstindex + 1)),
1956 gp_Pnt2d(LastParams.Value(afirstindex),
1957 LastParams.Value(afirstindex + 1)));
1960 while((d1.SquareMagnitude() < asqresol) &&
1961 (ii < aSeqOfNewPoint.Length()));
1965 while(nextii < aSeqOfNewPoint.Length()) {
1967 gp_Vec2d nextd1(0, 0);
1968 Standard_Integer jj = nextii;
1971 aSeqOfNewPoint(jj).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
1972 FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
1973 aSeqOfNewPoint(jj + 1).Parameters(LastParams.ChangeValue(0), LastParams.ChangeValue(1),
1974 LastParams.ChangeValue(2), LastParams.ChangeValue(3));
1975 nextd1 = gp_Vec2d(gp_Pnt2d(FirstParams.Value(afirstindex),
1976 FirstParams.Value(afirstindex + 1)),
1977 gp_Pnt2d(LastParams.Value(afirstindex),
1978 LastParams.Value(afirstindex + 1)));
1982 while((nextd1.SquareMagnitude() < asqresol) &&
1983 (jj < aSeqOfNewPoint.Length()));
1986 if(fabs(d1.Angle(nextd1)) > piquota) {
1987 bExtendLine = Standard_False;
1993 // end if(fabs(aTangentZoneDir.Angle(anIsoDir)
2000 return Standard_False;
2002 Standard_Integer i = 0;
2004 for(i = 1; i <= aSeqOfNewPoint.Length(); i++) {
2005 AddAPoint(line, aSeqOfNewPoint.Value(i));
2008 return bOutOfTangentZone;
2011 //=======================================================================
2012 //function : DistanceMinimizeByGradient
2016 // theInit should be initialized before function calling.
2017 //=======================================================================
2018 Standard_Boolean IntWalk_PWalking::
2019 DistanceMinimizeByGradient( const Handle(Adaptor3d_HSurface)& theASurf1,
2020 const Handle(Adaptor3d_HSurface)& theASurf2,
2021 TColStd_Array1OfReal& theInit,
2022 const Standard_Real* theStep0)
2024 const Standard_Integer aNbIterMAX = 60;
2025 const Standard_Real aTol = 1.0e-14;
2026 const Standard_Real aTolNul = 1.0 / Precision::Infinite();
2028 // I.e. if theU1 = 0.0 then Epsilon(theU1) = DBL_MIN (~1.0e-308).
2029 // Work with this number is impossible: there is a dangerous to
2030 // obtain Floating-point-overflow. Therefore, we limit this value.
2031 const Standard_Real aMinAddValU1 = Max(Epsilon(theInit(1)), aTolNul);
2032 const Standard_Real aMinAddValV1 = Max(Epsilon(theInit(2)), aTolNul);
2033 const Standard_Real aMinAddValU2 = Max(Epsilon(theInit(3)), aTolNul);
2034 const Standard_Real aMinAddValV2 = Max(Epsilon(theInit(4)), aTolNul);
2036 Handle(Geom_Surface) aS1, aS2;
2038 if (theASurf1->GetType() != GeomAbs_BezierSurface &&
2039 theASurf1->GetType() != GeomAbs_BSplineSurface)
2040 return Standard_True;
2041 if (theASurf2->GetType() != GeomAbs_BezierSurface &&
2042 theASurf2->GetType() != GeomAbs_BSplineSurface)
2043 return Standard_True;
2045 Standard_Boolean aStatus = Standard_False;
2048 gp_Vec aD1u, aD1v, aD2U, aD2V;
2050 theASurf1->D1(theInit(1), theInit(2), aP1, aD1u, aD1v);
2051 theASurf2->D1(theInit(3), theInit(4), aP2, aD2U, aD2V);
2053 Standard_Real aSQDistPrev = aP1.SquareDistance(aP2);
2055 gp_Vec aP12(aP1, aP2);
2057 Standard_Real aGradFu(-aP12.Dot(aD1u));
2058 Standard_Real aGradFv(-aP12.Dot(aD1v));
2059 Standard_Real aGradFU( aP12.Dot(aD2U));
2060 Standard_Real aGradFV( aP12.Dot(aD2V));
2062 Standard_Real aStepU1 = 1.0e-6, aStepV1 = 1.0e-6,
2063 aStepU2 = 1.0e-6, aStepV2 = 1.0e-6;
2067 aStepU1 = theStep0[0];
2068 aStepV1 = theStep0[1];
2069 aStepU2 = theStep0[2];
2070 aStepV2 = theStep0[3];
2073 Standard_Boolean flRepeat = Standard_True;
2074 Standard_Integer aNbIter = aNbIterMAX;
2078 Standard_Real anAdd = aGradFu*aStepU1;
2079 const Standard_Real aPARu = theInit(1) - Sign(Max(Abs(anAdd), aMinAddValU1), anAdd);
2081 anAdd = aGradFv*aStepV1;
2082 const Standard_Real aPARv = theInit(2) - Sign(Max(Abs(anAdd), aMinAddValV1), anAdd);
2084 anAdd = aGradFU*aStepU2;
2085 const Standard_Real aParU = theInit(3) - Sign(Max(Abs(anAdd), aMinAddValU2), anAdd);
2087 anAdd = aGradFV*aStepV2;
2088 const Standard_Real aParV = theInit(4) - Sign(Max(Abs(anAdd), aMinAddValV2), anAdd);
2092 theASurf1->D1(aPARu, aPARv, aPt1, aD1u, aD1v);
2093 theASurf2->D1(aParU, aParV, aPt2, aD2U, aD2V);
2095 Standard_Real aSQDist = aPt1.SquareDistance(aPt2);
2097 if(aSQDist < aSQDistPrev)
2099 aSQDistPrev = aSQDist;
2105 aStatus = aSQDistPrev < aTol;
2115 flRepeat = Standard_False;
2119 theASurf1->D1(theInit(1), theInit(2), aPt1, aD1u, aD1v);
2120 theASurf2->D1(theInit(3), theInit(4), aPt2, aD2U, aD2V);
2122 gp_Vec aPt12(aPt1, aPt2);
2123 aGradFu = -aPt12.Dot(aD1u);
2124 aGradFv = -aPt12.Dot(aD1v);
2125 aGradFU = aPt12.Dot(aD2U);
2126 aGradFV = aPt12.Dot(aD2V);
2130 aStepU1 = theStep0[0];
2131 aStepV1 = theStep0[1];
2132 aStepU2 = theStep0[2];
2133 aStepV2 = theStep0[3];
2137 aStepU1 = aStepV1 = aStepU2 = aStepV2 = 1.0e-6;
2146 //=======================================================================
2147 //function : DistanceMinimizeByExtrema
2151 // theP0, theU0 and theV0 parameters should be initialized
2152 // before the function calling.
2153 //=======================================================================
2154 Standard_Boolean IntWalk_PWalking::
2155 DistanceMinimizeByExtrema(const Handle(Adaptor3d_HSurface)& theASurf,
2156 const gp_Pnt& theP0,
2157 Standard_Real& theU0,
2158 Standard_Real& theV0,
2159 const Standard_Real* theStep0)
2161 const Standard_Real aTol = 1.0e-14;
2163 gp_Vec aD1Su, aD1Sv, aD2Su, aD2Sv, aD2SuvTemp;
2164 Standard_Real aSQDistPrev = RealLast();
2165 Standard_Real aU = theU0, aV = theV0;
2167 const Standard_Real aStep0[2] = { theStep0 ? theStep0[0] : 1.0,
2168 theStep0 ? theStep0[1] : 1.0 };
2170 Standard_Integer aNbIter = 10;
2173 theASurf->D2(aU, aV, aPS, aD1Su, aD1Sv, aD2Su, aD2Sv, aD2SuvTemp);
2175 gp_Vec aVec(theP0, aPS);
2177 Standard_Real aSQDist = aVec.SquareMagnitude();
2179 if(aSQDist >= aSQDistPrev)
2182 aSQDistPrev = aSQDist;
2187 if(aSQDistPrev < aTol)
2191 const Standard_Real aF1 = aD1Su.Dot(aVec), aF2 = aD1Sv.Dot(aVec);
2194 const Standard_Real aDf1u = aD2Su.Dot(aVec) + aD1Su.Dot(aD1Su),
2195 aDf1v = aD2Su.Dot(aD1Sv),
2197 aDf2v = aD2Sv.Dot(aVec) + aD1Sv.Dot(aD1Sv);
2199 const Standard_Real aDet = aDf1u*aDf2v - aDf1v*aDf2u;
2200 aU -= aStep0[0]*(aDf2v*aF1 - aDf1v*aF2) / aDet;
2201 aV += aStep0[1]*(aDf2u*aF1 - aDf1u*aF2) / aDet;
2205 return (aSQDistPrev < aTol);
2208 //=======================================================================
2209 //function : HandleSingleSingularPoint
2211 //=======================================================================
2212 Standard_Boolean IntWalk_PWalking::HandleSingleSingularPoint(const Handle(Adaptor3d_HSurface)& theASurf1,
2213 const Handle(Adaptor3d_HSurface)& theASurf2,
2214 const Standard_Real the3DTol,
2215 TColStd_Array1OfReal &thePnt)
2217 // u1, v1, u2, v2 order is used.
2218 Standard_Real aLowBorder[4] = {theASurf1->FirstUParameter(),
2219 theASurf1->FirstVParameter(),
2220 theASurf2->FirstUParameter(),
2221 theASurf2->FirstVParameter()};
2222 Standard_Real aUppBorder[4] = {theASurf1->LastUParameter(),
2223 theASurf1->LastVParameter(),
2224 theASurf2->LastUParameter(),
2225 theASurf2->LastVParameter()};
2226 IntImp_ConstIsoparametric aLockedDir[4] = {IntImp_UIsoparametricOnCaro1,
2227 IntImp_VIsoparametricOnCaro1,
2228 IntImp_UIsoparametricOnCaro2,
2229 IntImp_VIsoparametricOnCaro2};
2231 // Create new intersector with new tolerance.
2232 IntWalk_TheInt2S anInt(theASurf1, theASurf2, the3DTol);
2233 math_FunctionSetRoot aRsnld(anInt.Function());
2235 for (Standard_Integer i = 1; i <= 4; ++i)
2237 if ( Abs(thePnt(i) - aLowBorder[i - 1]) < Precision::PConfusion() ||
2238 Abs(thePnt(i) - aUppBorder[i - 1]) < Precision::PConfusion())
2241 anInt.Perform(thePnt,aRsnld, aLockedDir[i - 1]);
2243 if (!anInt.IsDone())
2246 if (anInt.IsEmpty())
2249 anInt.Point().Parameters(thePnt(1), thePnt(2), thePnt(3), thePnt(4));
2250 return Standard_True;
2254 return Standard_False;
2257 //=======================================================================
2258 //function : SeekPointOnBoundary
2260 //=======================================================================
2261 Standard_Boolean IntWalk_PWalking::
2262 SeekPointOnBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
2263 const Handle(Adaptor3d_HSurface)& theASurf2,
2264 const Standard_Real theU1,
2265 const Standard_Real theV1,
2266 const Standard_Real theU2,
2267 const Standard_Real theV2,
2268 const Standard_Boolean isTheFirst)
2270 Standard_Boolean isOK = Standard_False;
2272 // Tune solution tolerance according with object size.
2273 const Standard_Real aRes1 = Max(Precision::PConfusion() / theASurf1->UResolution(1.0),
2274 Precision::PConfusion() / theASurf1->VResolution(1.0));
2275 const Standard_Real aRes2 = Max(Precision::PConfusion() / theASurf2->UResolution(1.0),
2276 Precision::PConfusion() / theASurf2->VResolution(1.0));
2277 const Standard_Real a3DTol = Max(aRes1, aRes2);
2278 const Standard_Real aTol = Max(Precision::Confusion(), a3DTol);
2280 // u1, v1, u2, v2 order is used.
2281 TColStd_Array1OfReal aPnt(1,4);
2282 aPnt(1) = theU1; aPnt(2) = theV1; aPnt(3) = theU2; aPnt(4) = theV2;
2283 TColStd_Array1OfReal aSingularPnt(aPnt);
2285 Standard_Integer aNbIter = 20;
2286 Standard_Boolean aStatus = Standard_False;
2290 aStatus = DistanceMinimizeByGradient(theASurf1, theASurf2, aPnt);
2294 aStatus = DistanceMinimizeByExtrema(theASurf1, theASurf2->Value(aPnt(3), aPnt(4)), aPnt(1), aPnt(2));
2298 aStatus = DistanceMinimizeByExtrema(theASurf2, theASurf1->Value(aPnt(1), aPnt(2)), aPnt(3), aPnt(4));
2302 while(!aStatus && (aNbIter > 0));
2304 // Handle singular points.
2305 Standard_Boolean aSingularStatus = HandleSingleSingularPoint(theASurf1, theASurf2, aTol, aSingularPnt);
2306 if (aSingularStatus)
2307 aPnt = aSingularPnt;
2309 if (!aStatus && !aSingularStatus)
2314 gp_Pnt aP1 = theASurf1->Value(aPnt(1), aPnt(2));
2315 gp_Pnt aP2 = theASurf2->Value(aPnt(3), aPnt(4));
2316 const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2318 const Standard_Real aSQDist = aPInt.SquareDistance(aP1);
2319 if (aSQDist > aTol * aTol)
2324 //Found point is true intersection point
2325 IntSurf_PntOn2S anIP;
2326 anIP.SetValue(aPInt, aPnt(1), aPnt(2), aPnt(3), aPnt(4));
2328 //The main idea of checks below is to define if insertion of
2329 //addition point (on the boundary) does not lead to invalid
2330 //intersection curve (e.g. having a loop).
2332 //Loops are detected with rotation angle of the Walking-line (WL).
2333 //If there is hairpin bend then insertion is forbidden.
2335 //There are at least two possible problems:
2336 // 1. There are some cases when two neighbor points of the WL
2337 // are almost coincident (the distance between them is less
2338 // than Precision::Confusion). It is impossible to define
2339 // rotation angle in these cases. Therefore, points with
2340 // "good" distances should be selected.
2342 // 2. Intersection point on the surface boundary has highest
2343 // priority in compare with other "middle" points. Therefore,
2344 // if insertion of new point will result in a bend then some
2345 // "middle" points should be deleted in order to provide
2346 // correct insertion.
2348 //Problem test cases:
2349 // test bugs modalg_5 bug24585_1
2350 // test boolean bcut_complex G7
2351 // test bugs moddata_2 bug469
2355 while (line->NbPoints() > 1)
2357 const Standard_Integer aNbPnts = line->NbPoints();
2359 Standard_Integer aPInd = 1;
2360 for (; aPInd <= aNbPnts; aPInd++)
2362 aP1.SetXYZ(line->Value(aPInd).Value().XYZ());
2363 if (aP1.SquareDistance(aPInt) > Precision::SquareConfusion())
2367 for (++aPInd; aPInd <= aNbPnts; aPInd++)
2369 aP2.SetXYZ(line->Value(aPInd).Value().XYZ());
2370 if (aP1.SquareDistance(aP2) > Precision::SquareConfusion())
2374 if (aPInd > aNbPnts)
2379 const gp_XYZ aDir01(aP1.XYZ() - aPInt.XYZ());
2380 const gp_XYZ aDir12(aP2.XYZ() - aP1.XYZ());
2382 if (aDir01.Dot(aDir12) > 0.0)
2387 line->RemovePoint(1);
2390 line->InsertBefore(1, anIP);
2391 isOK = Standard_True;
2395 while (line->NbPoints() > 1)
2397 const Standard_Integer aNbPnts = line->NbPoints();
2399 gp_Pnt aPPrev, aPCurr;
2400 Standard_Integer aPInd = aNbPnts;
2401 for (; aPInd > 0; aPInd--)
2403 aPCurr.SetXYZ(line->Value(aPInd).Value().XYZ());
2404 if (aPCurr.SquareDistance(aPInt) > Precision::SquareConfusion())
2408 for (--aPInd; aPInd > 0; aPInd--)
2410 aPPrev.SetXYZ(line->Value(aPInd).Value().XYZ());
2411 if (aPCurr.SquareDistance(aPPrev) > Precision::SquareConfusion())
2420 const gp_XYZ aDirPC(aPCurr.XYZ() - aPPrev.XYZ());
2421 const gp_XYZ aDirCN(aPInt.XYZ() - aPCurr.XYZ());
2423 if (aDirPC.Dot(aDirCN) > 0.0)
2428 line->RemovePoint(aNbPnts);
2432 isOK = Standard_True;
2438 //=======================================================================
2439 //function : PutToBoundary
2441 //=======================================================================
2442 Standard_Boolean IntWalk_PWalking::
2443 PutToBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
2444 const Handle(Adaptor3d_HSurface)& theASurf2)
2446 const Standard_Real aTolMin = Precision::Confusion();
2448 Standard_Boolean hasBeenAdded = Standard_False;
2450 const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2451 const Standard_Real aU1bLast = theASurf1->LastUParameter();
2452 const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2453 const Standard_Real aU2bLast = theASurf2->LastUParameter();
2454 const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2455 const Standard_Real aV1bLast = theASurf1->LastVParameter();
2456 const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2457 const Standard_Real aV2bLast = theASurf2->LastVParameter();
2459 Standard_Real aTol = 1.0;
2460 aTol = Min(aTol, aU1bLast - aU1bFirst);
2461 aTol = Min(aTol, aU2bLast - aU2bFirst);
2462 aTol = Min(aTol, aV1bLast - aV1bFirst);
2463 aTol = Min(aTol, aV2bLast - aV2bFirst)*1.0e-3;
2465 if(aTol <= 2.0*aTolMin)
2466 return hasBeenAdded;
2468 Standard_Boolean isNeedAdding = Standard_False;
2469 Standard_Boolean isU1parallel = Standard_False, isV1parallel = Standard_False;
2470 Standard_Boolean isU2parallel = Standard_False, isV2parallel = Standard_False;
2471 IsParallel(line, Standard_True, aTol, isU1parallel, isV1parallel);
2472 IsParallel(line, Standard_False, aTol, isU2parallel, isV2parallel);
2474 Standard_Real u1, v1, u2, v2;
2475 line->Value(1).Parameters(u1, v1, u2, v2);
2476 Standard_Real aDelta = 0.0;
2480 aDelta = u1 - aU1bFirst;
2481 if((aTolMin < aDelta) && (aDelta < aTol))
2484 isNeedAdding = Standard_True;
2488 aDelta = aU1bLast - u1;
2489 if((aTolMin < aDelta) && (aDelta < aTol))
2492 isNeedAdding = Standard_True;
2499 aDelta = u2 - aU2bFirst;
2500 if((aTolMin < aDelta) && (aDelta < aTol))
2503 isNeedAdding = Standard_True;
2507 aDelta = aU2bLast - u2;
2508 if((aTolMin < aDelta) && (aDelta < aTol))
2511 isNeedAdding = Standard_True;
2518 aDelta = v1 - aV1bFirst;
2519 if((aTolMin < aDelta) && (aDelta < aTol))
2522 isNeedAdding = Standard_True;
2526 aDelta = aV1bLast - v1;
2527 if((aTolMin < aDelta) && (aDelta < aTol))
2530 isNeedAdding = Standard_True;
2537 aDelta = v2 - aV2bFirst;
2538 if((aTolMin < aDelta) && (aDelta < aTol))
2541 isNeedAdding = Standard_True;
2545 aDelta = aV2bLast - v2;
2546 if((aTolMin < aDelta) && (aDelta < aTol))
2549 isNeedAdding = Standard_True;
2557 SeekPointOnBoundary(theASurf1, theASurf2, u1,
2558 v1, u2, v2, Standard_True);
2561 const Standard_Integer aNbPnts = line->NbPoints();
2562 isNeedAdding = Standard_False;
2563 line->Value(aNbPnts).Parameters(u1, v1, u2, v2);
2567 aDelta = u1 - aU1bFirst;
2568 if((aTolMin < aDelta) && (aDelta < aTol))
2571 isNeedAdding = Standard_True;
2575 aDelta = aU1bLast - u1;
2576 if((aTolMin < aDelta) && (aDelta < aTol))
2579 isNeedAdding = Standard_True;
2586 aDelta = u2 - aU2bFirst;
2587 if((aTolMin < aDelta) && (aDelta < aTol))
2590 isNeedAdding = Standard_True;
2594 aDelta = aU2bLast - u2;
2595 if((aTolMin < aDelta) && (aDelta < aTol))
2598 isNeedAdding = Standard_True;
2605 aDelta = v1 - aV1bFirst;
2606 if((aTolMin < aDelta) && (aDelta < aTol))
2609 isNeedAdding = Standard_True;
2613 aDelta = aV1bLast - v1;
2614 if((aTolMin < aDelta) && (aDelta < aTol))
2617 isNeedAdding = Standard_True;
2624 aDelta = v2 - aV2bFirst;
2625 if((aTolMin < aDelta) && (aDelta < aTol))
2628 isNeedAdding = Standard_True;
2632 aDelta = aV2bLast - v2;
2633 if((aTolMin < aDelta) && (aDelta < aTol))
2636 isNeedAdding = Standard_True;
2644 SeekPointOnBoundary(theASurf1, theASurf2, u1,
2645 v1, u2, v2, Standard_False);
2648 return hasBeenAdded;
2651 //=======================================================================
2652 //function : SeekAdditionalPoints
2654 //=======================================================================
2655 Standard_Boolean IntWalk_PWalking::
2656 SeekAdditionalPoints( const Handle(Adaptor3d_HSurface)& theASurf1,
2657 const Handle(Adaptor3d_HSurface)& theASurf2,
2658 const Standard_Integer theMinNbPoints)
2660 const Standard_Real aTol = 1.0e-14;
2661 Standard_Integer aNbPoints = line->NbPoints();
2662 if(aNbPoints > theMinNbPoints)
2663 return Standard_True;
2665 const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2666 const Standard_Real aU1bLast = theASurf1->LastUParameter();
2667 const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2668 const Standard_Real aU2bLast = theASurf2->LastUParameter();
2669 const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2670 const Standard_Real aV1bLast = theASurf1->LastVParameter();
2671 const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2672 const Standard_Real aV2bLast = theASurf2->LastVParameter();
2675 Standard_Boolean isPrecise = Standard_False;
2677 TColStd_Array1OfReal aPnt(1, 4);
2680 Standard_Integer aNbPointsPrev = 0;
2681 while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
2683 aNbPointsPrev = aNbPoints;
2684 for(Standard_Integer fp = 1, lp = 2; fp < aNbPoints; fp = lp + 1)
2686 Standard_Real U1f, V1f, U2f, V2f; //first point in 1st and 2nd surafaces
2687 Standard_Real U1l, V1l, U2l, V2l; //last point in 1st and 2nd surafaces
2690 line->Value(fp).Parameters(U1f, V1f, U2f, V2f);
2691 line->Value(lp).Parameters(U1l, V1l, U2l, V2l);
2693 aPnt(1) = 0.5*(U1f + U1l);
2694 if(aPnt(1) < aU1bFirst)
2695 aPnt(1) = aU1bFirst;
2696 if(aPnt(1) > aU1bLast)
2699 aPnt(2) = 0.5*(V1f+V1l);
2700 if(aPnt(2) < aV1bFirst)
2701 aPnt(2) = aV1bFirst;
2702 if(aPnt(2) > aV1bLast)
2705 aPnt(3) = 0.5*(U2f+U2l);
2706 if(aPnt(3) < aU2bFirst)
2707 aPnt(3) = aU2bFirst;
2708 if(aPnt(3) > aU2bLast)
2711 aPnt(4) = 0.5*(V2f+V2l);
2712 if(aPnt(4) < aV2bFirst)
2713 aPnt(4) = aV2bFirst;
2714 if(aPnt(4) > aV2bLast)
2717 Standard_Boolean aStatus = Standard_False;
2718 Standard_Integer aNbIter = 5;
2721 aStatus = DistanceMinimizeByGradient(theASurf1, theASurf2, aPnt);
2727 aStatus = DistanceMinimizeByExtrema(theASurf1, theASurf2->Value(aPnt(3), aPnt(4)), aPnt(1), aPnt(2));
2733 aStatus = DistanceMinimizeByExtrema(theASurf2, theASurf1->Value(aPnt(1), aPnt(2)), aPnt(3), aPnt(4));
2739 while(!aStatus && (--aNbIter > 0));
2743 gp_Pnt aP1 = theASurf1->Value(aPnt(1), aPnt(2)),
2744 aP2 = theASurf2->Value(aPnt(3), aPnt(4));
2745 gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2747 const Standard_Real aSQDist1 = aPInt.SquareDistance(aP1),
2748 aSQDist2 = aPInt.SquareDistance(aP2);
2750 if((aSQDist1 < aTol) && (aSQDist2 < aTol))
2752 IntSurf_PntOn2S anIP;
2753 anIP.SetValue(aPInt, aPnt(1), aPnt(2), aPnt(3), aPnt(4));
2754 line->InsertBefore(lp, anIP);
2756 isPrecise = Standard_True;
2758 if(++aNbPoints >= theMinNbPoints)
2772 void IntWalk_PWalking::
2773 RepartirOuDiviser(Standard_Boolean& DejaReparti,
2774 IntImp_ConstIsoparametric& ChoixIso,
2775 Standard_Boolean& Arrive)
2777 // at the neighborhood of a point, there is a fail of marching
2778 // it is required to divide the steps to try to continue
2779 // if the step is too small if we are on border
2780 // restart in another direction if it was not done, otherwise stop
2783 // Standard_Integer i;
2784 if (Arrive) { //restart in the other direction
2785 if (!DejaReparti ) {
2786 Arrive = Standard_False;
2787 DejaReparti = Standard_True;
2788 previousPoint = line->Value(1);
2789 previoustg = Standard_False;
2790 previousd1 = firstd1;
2791 previousd2 = firstd2;
2793 indextg = line->NbPoints();
2797 //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2798 sensCheminement = -1;
2800 tglast = Standard_False;
2801 ChoixIso = choixIsoSav;
2808 Standard_Real u1,v1,u2,v2;
2809 Standard_Real U1,V1,U2,V2;
2810 Standard_Integer nn=line->NbPoints();
2812 line->Value(nn).Parameters(u1,v1,u2,v2);
2813 line->Value(nn-1).Parameters(U1,V1,U2,V2);
2814 pasuv[0]=Abs(u1-U1);
2815 pasuv[1]=Abs(v1-V1);
2816 pasuv[2]=Abs(u2-U2);
2817 pasuv[3]=Abs(v2-V2);
2824 if ( pasuv[0]*0.5 < ResoU1
2825 && pasuv[1]*0.5 < ResoV1
2826 && pasuv[2]*0.5 < ResoU2
2827 && pasuv[3]*0.5 < ResoV2
2830 tglast = Standard_True; // IS IT ENOUGH ????
2833 if (!DejaReparti) { //restart in the other direction
2834 DejaReparti = Standard_True;
2835 previousPoint = line->Value(1);
2836 previoustg = Standard_False;
2837 previousd1 = firstd1;
2838 previousd2 = firstd2;
2840 indextg = line->NbPoints();
2844 //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2846 sensCheminement = -1;
2848 tglast = Standard_False;
2849 ChoixIso = choixIsoSav;
2857 Standard_Real u1,v1,u2,v2;
2858 Standard_Real U1,V1,U2,V2;
2859 Standard_Integer nn=line->NbPoints();
2861 line->Value(nn).Parameters(u1,v1,u2,v2);
2862 line->Value(nn-1).Parameters(U1,V1,U2,V2);
2863 pasuv[0]=Abs(u1-U1);
2864 pasuv[1]=Abs(v1-V1);
2865 pasuv[2]=Abs(u2-U2);
2866 pasuv[3]=Abs(v2-V2);
2870 else Arrive = Standard_True;
2882 //OCC431(apo): modified ->
2883 static const Standard_Real CosRef2D = Cos(M_PI/9.0), AngRef2D = M_PI/2.0;
2885 static const Standard_Real d = 7.0;
2888 IntWalk_StatusDeflection IntWalk_PWalking::TestDeflection(const IntImp_ConstIsoparametric choixIso)
2890 // test if vector is observed by calculating an increase of vector
2891 // or the previous point and its tangent, the new calculated point and its
2892 // tangent; it is possible to find a cube passing by the 2 points and having as a
2893 // derivative the tangents of the intersection
2894 // calculate the point with parameter 0.5 on cube=p1
2895 // calculate the medium point of 2 points of intersection=p2
2896 // if arrow/2<=||p1p2||<= arrow consider that the vector is observed
2897 // otherwise adjust the step depending on the ratio ||p1p2||/vector
2898 // and the previous step
2899 // test if in 2 tangent planes of surfaces there is no too great angle2d
2900 // grand : if yes divide the step
2901 // test if there is no change of side
2904 if(line->NbPoints() ==1 ) {
2905 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=STATIC_PRECEDENT_INFLEXION=0;
2908 IntWalk_StatusDeflection Status = IntWalk_OK;
2909 Standard_Real FlecheCourante , Ratio = 1.0;
2912 const Handle(Adaptor3d_HSurface)& Caro1 = myIntersectionOn2S.Function().AuxillarSurface1();
2913 const Handle(Adaptor3d_HSurface)& Caro2 = myIntersectionOn2S.Function().AuxillarSurface2();
2915 const IntSurf_PntOn2S& CurrentPoint = myIntersectionOn2S.Point();
2916 //==================================================================================
2917 //========= S t o p o n p o i n t ============
2918 //==================================================================================
2919 if (myIntersectionOn2S.IsTangent()) {
2920 return IntWalk_ArretSurPoint;
2923 const gp_Dir& TgCourante = myIntersectionOn2S.Direction();
2925 const Standard_Real aCosBetweenTangent = TgCourante.Dot(previousd);
2927 //==================================================================================
2928 //========= R i s k o f i n f l e x i o n p o i n t ============
2929 //==================================================================================
2930 if (aCosBetweenTangent < 0) {
2931 //------------------------------------------------------------
2932 //-- Risk of inflexion point : Divide the step by 2
2933 //-- Initialize STATIC_PRECEDENT_INFLEXION so that
2934 //-- at the next call to return Pas_OK if there is no
2935 //-- more risk of the point of inflexion
2936 //------------------------------------------------------------
2942 STATIC_PRECEDENT_INFLEXION+=3;
2943 if (pasuv[0] < ResoU1 && pasuv[1] <ResoV1 && pasuv[2] <ResoU2 && pasuv[3] < ResoV2)
2944 return IntWalk_ArretSurPointPrecedent;
2946 return IntWalk_PasTropGrand;
2949 if(STATIC_PRECEDENT_INFLEXION > 0) {
2950 STATIC_PRECEDENT_INFLEXION -- ;
2955 //==================================================================================
2956 //========= D e t e c t c o n f u s e d P o in t s ===========
2957 //==================================================================================
2959 const Standard_Real aSqDist = previousPoint.Value().
2960 SquareDistance(CurrentPoint.Value());
2963 if (aSqDist < Precision::SquareConfusion()) {
2964 pasInit[0] = Max(pasInit[0], 5.0*ResoU1);
2965 pasInit[1] = Max(pasInit[1], 5.0*ResoV1);
2966 pasInit[2] = Max(pasInit[2], 5.0*ResoU2);
2967 pasInit[3] = Max(pasInit[3], 5.0*ResoV2);
2969 for(Standard_Integer i = 0; i < 4; i++)
2971 pasuv[i] = Max(pasuv[i], Min(1.5*pasuv[i], pasInit[i]));
2973 //Compute local resolution: for OCC26717
2974 if (Abs(pasuv[choixIso] - pasInit[choixIso]) <= Precision::Confusion())
2976 Standard_Real CurU, CurV;
2977 if (choixIso == IntImp_UIsoparametricOnCaro1 ||
2978 choixIso == IntImp_VIsoparametricOnCaro1)
2979 previousPoint.ParametersOnS1(CurU, CurV);
2981 previousPoint.ParametersOnS2(CurU, CurV);
2982 gp_Pnt CurPnt = (choixIso == IntImp_UIsoparametricOnCaro1 ||
2983 choixIso == IntImp_VIsoparametricOnCaro1)?
2984 Adaptor3d_HSurfaceTool::Value(Caro1, CurU, CurV) :
2985 Adaptor3d_HSurfaceTool::Value(Caro2, CurU, CurV);
2989 case IntImp_UIsoparametricOnCaro1:
2991 Adaptor3d_HSurfaceTool::Value(Caro1,
2992 CurU + sensCheminement*pasuv[0],
2995 case IntImp_VIsoparametricOnCaro1:
2997 Adaptor3d_HSurfaceTool::Value(Caro1,
2999 CurV + sensCheminement*pasuv[1]);
3001 case IntImp_UIsoparametricOnCaro2:
3003 Adaptor3d_HSurfaceTool::Value(Caro2,
3004 CurU + sensCheminement*pasuv[2],
3007 case IntImp_VIsoparametricOnCaro2:
3009 Adaptor3d_HSurfaceTool::Value(Caro2,
3011 CurV + sensCheminement*pasuv[3]);
3015 Standard_Real RefDist = CurPnt.Distance(OffsetPnt);
3016 Standard_Real LocalResol = 0.;
3017 if (RefDist > gp::Resolution())
3018 LocalResol = pasuv[choixIso] * tolconf / RefDist;
3019 if (pasuv[choixIso] < 2*LocalResol)
3020 pasuv[choixIso] = pasInit[choixIso] = 2*LocalResol;
3022 ////////////////////////////////////////
3023 Status = IntWalk_PointConfondu;
3026 //==================================================================================
3027 Standard_Real Up1,Vp1,Uc1,Vc1,Du1,Dv1,AbsDu1,AbsDu2,AbsDv1,AbsDv2;
3028 Standard_Real Up2,Vp2,Uc2,Vc2,Du2,Dv2;
3030 previousPoint.Parameters(Up1,Vp1,Up2,Vp2);
3031 CurrentPoint.Parameters(Uc1,Vc1,Uc2,Vc2);
3033 Du1 = Uc1 - Up1; Dv1 = Vc1 - Vp1;
3034 Du2 = Uc2 - Up2; Dv2 = Vc2 - Vp2;
3040 //=================================================================================
3041 //==== S t e p o f p r o g r e s s i o n (between previous and Current) =======
3042 //=================================================================================
3043 if ( AbsDu1 < ResoU1 && AbsDv1 < ResoV1
3044 && AbsDu2 < ResoU2 && AbsDv2 < ResoV2) {
3045 pasuv[0] = ResoU1; pasuv[1] = ResoV1; pasuv[2] = ResoU2; pasuv[3] = ResoV2;
3046 return(IntWalk_ArretSurPointPrecedent);
3048 //==================================================================================
3050 Standard_Real tolArea = 100.0;
3051 if (ResoU1 < Precision::PConfusion() ||
3052 ResoV1 < Precision::PConfusion() ||
3053 ResoU2 < Precision::PConfusion() ||
3054 ResoV2 < Precision::PConfusion() )
3055 tolArea = tolArea*2.0;
3057 Standard_Real Cosi1, CosRef1, Ang1, AngRef1, ResoUV1, Duv1, d1, tolCoeff1;
3058 Standard_Real Cosi2, CosRef2, Ang2, AngRef2, ResoUV2, Duv2, d2, tolCoeff2;
3059 Cosi1 = Du1*previousd1.X() + Dv1*previousd1.Y();
3060 Cosi2 = Du2*previousd2.X() + Dv2*previousd2.Y();
3061 Duv1 = Du1*Du1 + Dv1*Dv1;
3062 Duv2 = Du2*Du2 + Dv2*Dv2;
3063 ResoUV1 = ResoU1*ResoU1 + ResoV1*ResoV1;
3064 ResoUV2 = ResoU2*ResoU2 + ResoV2*ResoV2;
3066 //modified by NIZNHY-PKV Wed Nov 13 12:25:44 2002 f
3068 Standard_Real aMinDiv2=Precision::Confusion();
3069 aMinDiv2=aMinDiv2*aMinDiv2;
3072 if (Duv1>aMinDiv2) {
3073 d1 = Abs(ResoUV1/Duv1);
3074 d1 = Min(Sqrt(d1)*tolArea, d);
3076 //d1 = Abs(ResoUV1/Duv1);
3077 //d1 = Min(Sqrt(d1)*tolArea,d);
3078 //modified by NIZNHY-PKV Wed Nov 13 12:34:30 2002 t
3079 tolCoeff1 = Exp(d1);
3081 //modified by NIZNHY-PKV Wed Nov 13 12:34:43 2002 f
3083 if (Duv2>aMinDiv2) {
3084 d2 = Abs(ResoUV2/Duv2);
3085 d2 = Min(Sqrt(d2)*tolArea,d);
3087 //d2 = Abs(ResoUV2/Duv2);
3088 //d2 = Min(Sqrt(d2)*tolArea,d);
3089 //modified by NIZNHY-PKV Wed Nov 13 12:34:53 2002 t
3090 tolCoeff2 = Exp(d2);
3091 CosRef1 = CosRef2D/tolCoeff1;
3092 CosRef2 = CosRef2D/tolCoeff2;
3094 //==================================================================================
3095 //== The points are not confused : ==
3096 //== D e t e c t t h e S t o p a t p r e v i o u s p o i n t ==
3097 //== N o t T o o G r e a t (angle in space UV) ==
3098 //== C h a n g e o f s i d e ==
3099 //==================================================================================
3100 if (Status != IntWalk_PointConfondu) {
3101 if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2) {
3102 pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3103 if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2) {
3104 return(IntWalk_ArretSurPointPrecedent);
3107 pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3108 return(IntWalk_PasTropGrand);
3111 const gp_Dir2d& Tg2dcourante1 = myIntersectionOn2S.DirectionOnS1();
3112 const gp_Dir2d& Tg2dcourante2 = myIntersectionOn2S.DirectionOnS2();
3113 Cosi1 = Du1*Tg2dcourante1.X() + Dv1*Tg2dcourante1.Y();
3114 Cosi2 = Du2*Tg2dcourante2.X() + Dv2*Tg2dcourante2.Y();
3115 Ang1 = Abs(previousd1.Angle(Tg2dcourante1));
3116 Ang2 = Abs(previousd2.Angle(Tg2dcourante2));
3117 AngRef1 = AngRef2D*tolCoeff1;
3118 AngRef2 = AngRef2D*tolCoeff2;
3119 //-------------------------------------------------------
3120 //-- Test : Angle too great in space UV -----
3121 //-- Change of side -----
3122 //-------------------------------------------------------
3123 if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2 || Ang1 > AngRef1 || Ang2 > AngRef2) {
3124 pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3125 if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2)
3126 return(IntWalk_ArretSurPoint);
3128 return(IntWalk_PasTropGrand);
3132 //==================================================================================
3133 //== D e t e c t i o n o f : Step Too Small
3135 //==================================================================================
3137 //---------------------------------------
3138 //-- Estimate of the vector --
3139 //---------------------------------------
3141 Sqrt(Abs((previousd.XYZ()-TgCourante.XYZ()).SquareModulus()*aSqDist))/8.;
3143 if ( FlecheCourante<= fleche*0.5) { //-- Current step too small
3144 if(FlecheCourante>1e-16) {
3145 Ratio = 0.5*(fleche/FlecheCourante);
3150 Standard_Real pasSu1 = pasuv[0];
3151 Standard_Real pasSv1 = pasuv[1];
3152 Standard_Real pasSu2 = pasuv[2];
3153 Standard_Real pasSv2 = pasuv[3];
3156 //-- a point at U+DeltaU is required, ....
3157 //-- return a point at U + Epsilon
3158 //-- Epsilon << DeltaU.
3160 if(pasuv[0]< AbsDu1) pasuv[0] = AbsDu1;
3161 if(pasuv[1]< AbsDv1) pasuv[1] = AbsDv1;
3162 if(pasuv[2]< AbsDu2) pasuv[2] = AbsDu2;
3163 if(pasuv[3]< AbsDv2) pasuv[3] = AbsDv2;
3165 if(pasuv[0]<ResoU1) pasuv[0]=ResoU1;
3166 if(pasuv[1]<ResoV1) pasuv[1]=ResoV1;
3167 if(pasuv[2]<ResoU2) pasuv[2]=ResoU2;
3168 if(pasuv[3]<ResoV2) pasuv[3]=ResoV2;
3169 //-- if(Ratio>10.0 ) { Ratio=10.0; }
3170 Standard_Real R1,R = pasInit[0]/pasuv[0];
3171 R1= pasInit[1]/pasuv[1]; if(R1<R) R=R1;
3172 R1= pasInit[2]/pasuv[2]; if(R1<R) R=R1;
3173 R1= pasInit[3]/pasuv[3]; if(R1<R) R=R1;
3174 if(Ratio > R) Ratio=R;
3175 pasuv[0] = Min(Ratio*pasuv[0],pasInit[0]);
3176 pasuv[1] = Min(Ratio*pasuv[1],pasInit[1]);
3177 pasuv[2] = Min(Ratio*pasuv[2],pasInit[2]);
3178 pasuv[3] = Min(Ratio*pasuv[3],pasInit[3]);
3179 if (pasuv[0] != pasSu1 || pasuv[2] != pasSu2||
3180 pasuv[1] != pasSv1 || pasuv[3] != pasSv2) {
3181 if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3182 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3183 return IntWalk_PasTropGrand;
3186 if(Status == IntWalk_OK) {
3187 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
3188 //-- Try to increase the step
3192 else { //-- CurrentVector > vector*0.5
3193 if (FlecheCourante > fleche) { //-- Current step too Great
3194 Ratio = fleche/FlecheCourante;
3195 pasuv[0] = Ratio*pasuv[0];
3196 pasuv[1] = Ratio*pasuv[1];
3197 pasuv[2] = Ratio*pasuv[2];
3198 pasuv[3] = Ratio*pasuv[3];
3199 //if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3200 // STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3201 return IntWalk_PasTropGrand;
3204 else { //-- vector/2 < CurrentVector <= vector
3205 Ratio = 0.75 * (fleche / FlecheCourante);
3209 if(Status != IntWalk_PointConfondu)
3211 //Here, aCosBetweenTangent >= 0.0 definitely.
3214 Brief algorithm description.
3215 We have two (not-coincindent) intersection point (P1 and P2). In every point,
3216 vector of tangent (to the intersection curve) is known (vectors T1 and T2).
3217 Basing on these data, we create osculating circle.
3219 * - arc of osculating circle
3226 Let me draw your attention to the following facts:
3227 1. Vectors T1 and T2 direct FROM (not TO) points P1 and P2. Therefore,
3228 one of previously computed vector should be reversed.
3230 In this case, the absolute (!) value of the deflection between the arc of
3231 the osculating circle and the P1P2 segment can be computed as follows:
3232 e = d*(1-sin(B/2))/(2*cos(B/2)), (1)
3233 where d is the length of P1P2 segment, B is the angle between vectors T1 and T2.
3241 Later, the normal state of algorithm work is (as we apply)
3242 tolconf/2 <= e <= tolconf.
3243 In this case, we keep previous step.
3245 If e < tolconf/2 then the local curvature of the intersection curve is small.
3246 As result, the step should be increased.
3248 If e > tolconf then the step is too big. Therefore, we should decrease one.
3250 Condition (1) is equivalent to
3251 sin(B/2) = 1 - 2/(1+(d/(2*e))^2) = Fs(e),
3252 cos(B) = 1 - 2*Fs(e)^2 = Fd(e),
3253 where Fs(e)and Fd(e) are some function with parameter "deflection".
3255 Let mean that Fs(e) is decreasing function. Fd(e) is increasing function,
3256 in the range, where Fs(e) > 0.0 (i.e. when e < d/2).
3258 Now, let substitute required deflection (tolconf or tolconf/2) to z. Then
3259 it is necessary to check if e < z or if e > z.
3261 In this case, it is enough to comapare Fs(e) and Fs(z).
3262 At that Fs(e) > 0 because sin(B/2) > 0 always.
3264 Therefore, if Fs(z) < 0.0 then Fs(e) > Fs(z) ==> e < z definitely.
3265 If Fs(z) > 0.0 then we can compare Fs(z)^2 and Fs(e)^2 or, in substance,
3266 values Fd(e) and Fd(z). If Fd(e) > Fd(z) then e > z and vice versa.
3269 //Fd(e) is already known (Fd(e) == -aCosBetweenTangent)
3271 const Standard_Real anInvSqAbsArcDeflMax = 0.25*aSqDist/(tolconf*tolconf);
3272 const Standard_Real aSinB2Max = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMax);
3274 if(aSinB2Max >= 0.0 && (aCosBetweenTangent <= 2.0 * aSinB2Max * aSinB2Max - 1.0))
3275 {//Real deflection is greater or equal than tolconf
3276 Status = IntWalk_PasTropGrand;
3279 {//Real deflection is less than tolconf
3280 const Standard_Real anInvSqAbsArcDeflMin = 4.0*anInvSqAbsArcDeflMax;
3281 const Standard_Real aSinB2Min = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMin);
3283 if((aSinB2Min < 0.0) || (aCosBetweenTangent >= 2.0 * aSinB2Min * aSinB2Min - 1.0))
3284 {//Real deflection is less than tolconf/2.0
3285 Status = IntWalk_StepTooSmall;
3289 if(Status == IntWalk_PasTropGrand)
3291 pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3295 if(Status == IntWalk_StepTooSmall)
3297 pasuv[0] = Max(pasuv[0], AbsDu1);
3298 pasuv[1] = Max(pasuv[1], AbsDv1);
3299 pasuv[2] = Max(pasuv[2], AbsDu2);
3300 pasuv[3] = Max(pasuv[3], AbsDv2);
3302 pasInit[0] = Max(pasInit[0], AbsDu1);
3303 pasInit[1] = Max(pasInit[1], AbsDv1);
3304 pasInit[2] = Max(pasInit[2], AbsDu2);
3305 pasInit[3] = Max(pasInit[3], AbsDv2);
3311 pasuv[0] = Max(myStepMin[0],Min(Min(Ratio*AbsDu1,pasuv[0]),pasInit[0]));
3312 pasuv[1] = Max(myStepMin[1],Min(Min(Ratio*AbsDv1,pasuv[1]),pasInit[1]));
3313 pasuv[2] = Max(myStepMin[2],Min(Min(Ratio*AbsDu2,pasuv[2]),pasInit[2]));
3314 pasuv[3] = Max(myStepMin[3],Min(Min(Ratio*AbsDv2,pasuv[3]),pasInit[3]));
3316 if(Status == IntWalk_OK) STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
3320 Standard_Boolean IntWalk_PWalking::
3321 TestArret(const Standard_Boolean DejaReparti,
3322 TColStd_Array1OfReal& Param,
3323 IntImp_ConstIsoparametric& ChoixIso)
3326 // test if the point of intersection set by these parameters remains in the
3327 // natural domain of each square.
3328 // if the point outpasses reframe to find the best iso (border)
3329 // that intersects easiest the other square
3330 // otherwise test if closed line is present
3333 Standard_Real Uvd[4],Uvf[4],Epsuv[4],Duv[4],Uvp[4],dv,dv2,ParC[4];
3334 Standard_Real DPc,DPb;
3335 Standard_Integer i = 0, k = 0;
3340 previousPoint.Parameters(Uvp[0],Uvp[1],Uvp[2],Uvp[3]);
3342 Standard_Real SolParam[4];
3343 myIntersectionOn2S.Point().Parameters(SolParam[0],SolParam[1],SolParam[2],SolParam[3]);
3345 Standard_Boolean Trouve = Standard_False;
3347 Uvd[0]=Um1; Uvf[0]=UM1; Uvd[1]=Vm1; Uvf[1]=VM1;
3348 Uvd[2]=Um2; Uvf[2]=UM2; Uvd[3]=Vm2; Uvf[3]=VM2;
3350 Standard_Integer im1;
3351 for ( i = 1,im1 = 0;i<=4;i++,im1++) {
3358 if (Param(i) < (Uvd[im1]-Epsuv[im1]) ||
3359 SolParam[im1] < (Uvd[im1]-Epsuv[im1])) //-- Current ----- Bound Inf ----- Previous
3361 Trouve = Standard_True; //--
3362 DPc = Uvp[im1]-Param(i); //-- Previous - Current
3363 DPb = Uvp[im1]-Uvd[im1]; //-- Previous - Bound Inf
3364 ParC[im1] = Uvd[im1]; //-- ParamCorrige
3365 dv = Param(k)-Uvp[k-1]; //-- Current - Previous (other Direction)
3367 if(dv2>RealEpsilon()) { //-- Progress at the other Direction ?
3368 Duv[im1] = DPc*DPb + dv2;
3369 Duv[im1] = Duv[im1]*Duv[im1]/(DPc*DPc+dv2)/(DPb*DPb+dv2);
3372 Duv[im1]=-1.0; //-- If no progress, do not change
3373 } //-- the choice of iso
3375 else if (Param(i) > (Uvf[im1] + Epsuv[im1]) ||
3376 SolParam[im1] > (Uvf[im1] + Epsuv[im1]))//-- Previous ----- Bound Sup ----- Current
3378 Trouve = Standard_True; //--
3379 DPc = Param(i)-Uvp[im1]; //-- Current - Previous
3380 DPb = Uvf[im1]-Uvp[im1]; //-- Bound Sup - Previous
3381 ParC[im1] = Uvf[im1]; //-- Param Corrige
3382 dv = Param(k)-Uvp[k-1]; //-- Current - Previous (other Direction)
3384 if(dv2>RealEpsilon()) { //-- Progress in other Direction ?
3385 Duv[im1] = DPc*DPb + dv2;
3386 Duv[im1] = Duv[im1]*Duv[im1]/(DPc*DPc+dv2)/(DPb*DPb+dv2);
3389 Duv[im1]=-1.0; //-- If no progress, do not change
3390 } //-- the choice of iso
3399 //--------------------------------------------------
3400 //-- One of Parameters u1,v1,u2,v2 is outside of --
3401 //-- the natural limits. --
3402 //-- Find the best direction of --
3403 //-- progress and reframe the parameters. --
3404 //--------------------------------------------------
3405 Standard_Real ddv = -1.0;
3407 for (i=0;i<=3;i++) {
3408 Param(i+1) = ParC[i];
3415 ChoixIso = ChoixRef[k];
3418 if((ParC[0]<=Uvd[0]+Epsuv[0]) || (ParC[0]>=Uvf[0]-Epsuv[0])) {
3419 ChoixIso = IntImp_UIsoparametricOnCaro1;
3421 else if((ParC[1]<=Uvd[1]+Epsuv[1]) || (ParC[1]>=Uvf[1]-Epsuv[1])) {
3422 ChoixIso = IntImp_VIsoparametricOnCaro1;
3424 else if((ParC[2]<=Uvd[2]+Epsuv[2]) || (ParC[2]>=Uvf[2]-Epsuv[2])) {
3425 ChoixIso = IntImp_UIsoparametricOnCaro2;
3427 else if((ParC[3]<=Uvd[3]+Epsuv[3]) || (ParC[3]>=Uvf[3]-Epsuv[3])) {
3428 ChoixIso = IntImp_VIsoparametricOnCaro2;
3431 close = Standard_False;
3432 return Standard_True;
3436 if (!DejaReparti) { // find if line closed
3439 const IntSurf_PntOn2S& POn2S1=line->Value(1);
3441 POn2S1.ParametersOnS1(u,v);
3442 gp_Pnt2d P1uvS1(u,v);
3443 previousPoint.ParametersOnS1(u,v);
3444 gp_Pnt2d PrevuvS1(u,v);
3445 myIntersectionOn2S.Point().ParametersOnS1(u,v);
3446 gp_Pnt2d myIntersuvS1(u,v);
3447 Standard_Boolean close2dS1 = (P1uvS1.XY()-PrevuvS1.XY())*
3448 (P1uvS1.XY()-myIntersuvS1.XY()) < 0.0;
3450 POn2S1.ParametersOnS2(u,v);
3451 gp_Pnt2d P1uvS2(u,v);
3452 previousPoint.ParametersOnS2(u,v);
3453 gp_Pnt2d PrevuvS2(u,v);
3454 myIntersectionOn2S.Point().ParametersOnS2(u,v);
3455 gp_Pnt2d myIntersuvS2(u,v);
3456 Standard_Boolean close2dS2 = (P1uvS2.XY()-PrevuvS2.XY())*
3457 (P1uvS2.XY()-myIntersuvS2.XY()) < 0.0;
3459 close = close2dS1 && close2dS2;
3462 else return Standard_False;