Warning, /include/opencascade/IntStart_SearchOnBoundaries.gxx is written in an unsupported language. File is not indexed.
0001 // Copyright (c) 1995-1999 Matra Datavision
0002 // Copyright (c) 1999-2014 OPEN CASCADE SAS
0003 //
0004 // This file is part of Open CASCADE Technology software library.
0005 //
0006 // This library is free software; you can redistribute it and/or modify it under
0007 // the terms of the GNU Lesser General Public License version 2.1 as published
0008 // by the Free Software Foundation, with special exception defined in the file
0009 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0010 // distribution for complete text of the license and disclaimer of any warranty.
0011 //
0012 // Alternatively, this file may be used under the terms of Open CASCADE
0013 // commercial license or contractual agreement.
0014
0015 #include <algorithm>
0016 #include <memory>
0017 #include <TopoDS_Edge.hxx>
0018 #include <Geom_Curve.hxx>
0019 #include <BRepAdaptor_Curve.hxx>
0020 #include <Adaptor3d_Surface.hxx>
0021 #include <Adaptor3d_CurveOnSurface.hxx>
0022 #include <Adaptor3d_CurveOnSurface.hxx>
0023 #include <GeomAbs_SurfaceType.hxx>
0024 #include <BRep_Tool.hxx>
0025 #include <Geom_Line.hxx>
0026 #include <Geom_Plane.hxx>
0027 #include <Geom_CylindricalSurface.hxx>
0028 #include <Geom_ConicalSurface.hxx>
0029 #include <Geom_SphericalSurface.hxx>
0030 #include <Geom_ToroidalSurface.hxx>
0031 #include <gp_Lin.hxx>
0032 #include <gp_Vec.hxx>
0033 #include <gp_Dir.hxx>
0034 #include <gp_Cylinder.hxx>
0035 #include <gp_Ax1.hxx>
0036 #include <gp_Lin.hxx>
0037
0038 #include <GeomAdaptor_Curve.hxx>
0039 #include <GeomAdaptor_Surface.hxx>
0040 #include <Precision.hxx>
0041 #include <Extrema_ExtCC.hxx>
0042 //#include <Extrema_ExtCS.hxx>
0043 #include <Extrema_POnCurv.hxx>
0044 #include <IntCurveSurface_HInter.hxx>
0045
0046 #include <math_FunctionSample.hxx>
0047 #include <math_FunctionAllRoots.hxx>
0048 #include <TColgp_SequenceOfPnt.hxx>
0049
0050 // Modified by skv - Tue Aug 31 12:13:51 2004 OCC569
0051
0052 #include <Precision.hxx>
0053 #include <IntSurf_Quadric.hxx>
0054 #include <math_Function.hxx>
0055 #include <math_BrentMinimum.hxx>
0056 #include <math_Matrix.hxx>
0057 #include <math_Vector.hxx>
0058 #include <NCollection_Array1.hxx>
0059
0060 #ifdef OCCT_DEBUG
0061 #include <Geom_Circle.hxx>
0062 #include <Geom_Ellipse.hxx>
0063 #include <Geom_Hyperbola.hxx>
0064 #include <Geom_Parabola.hxx>
0065 #include <Geom_BezierCurve.hxx>
0066 #include <Geom_BSplineCurve.hxx>
0067 #include <GeomLib.hxx>
0068 #endif
0069
0070
0071 static Standard_Boolean IsDegenerated(const Handle(Adaptor3d_CurveOnSurface)& theCurve);
0072 static Standard_Boolean IsDegenerated(const IntSurf_Quadric& theQuadric);
0073
0074 static void FindVertex (const TheArc&,
0075 const Handle(TheTopolTool)&,
0076 TheFunction&,
0077 IntStart_SequenceOfPathPoint&,
0078 const Standard_Real);
0079
0080
0081 static void BoundedArc (const TheArc& A,
0082 const Handle(TheTopolTool)& Domain,
0083 const Standard_Real Pdeb,
0084 const Standard_Real Pfin,
0085 TheFunction& Func,
0086 IntStart_SequenceOfPathPoint& pnt,
0087 IntStart_SequenceOfSegment& seg,
0088 const Standard_Real TolBoundary,
0089 const Standard_Real TolTangency,
0090 Standard_Boolean& Arcsol,
0091 const Standard_Boolean RecheckOnRegularity);
0092
0093 static void PointProcess (const gp_Pnt&,
0094 const Standard_Real,
0095 const TheArc&,
0096 const Handle(TheTopolTool)&,
0097 IntStart_SequenceOfPathPoint&,
0098 const Standard_Real,
0099 Standard_Integer&);
0100
0101 static Standard_Integer TreatLC (const TheArc& A,
0102 const Handle(TheTopolTool)& aDomain,
0103 const IntSurf_Quadric& aQuadric,
0104 const Standard_Real TolBoundary,
0105 IntStart_SequenceOfPathPoint& pnt);
0106
0107 static Standard_Boolean IsRegularity(const TheArc& A,
0108 const Handle(TheTopolTool)& aDomain);
0109
0110 class MinFunction : public math_Function
0111 {
0112 public:
0113 MinFunction(TheFunction &theFunc) : myFunc(&theFunc) {};
0114
0115 //returns value of the one-dimension-function when parameter
0116 //is equal to theX
0117 virtual Standard_Boolean Value(const Standard_Real theX,
0118 Standard_Real& theFVal)
0119 {
0120 if(!myFunc->Value(theX, theFVal))
0121 return Standard_False;
0122
0123 theFVal *= theFVal;
0124 return Standard_True;
0125 }
0126
0127 //see analogical method for abstract owner class math_Function
0128 virtual Standard_Integer GetStateNumber()
0129 {
0130 return 0;
0131 }
0132
0133 private:
0134 TheFunction *myFunc;
0135 };
0136
0137
0138 //=======================================================================
0139 //function : FindVertex
0140 //purpose :
0141 //=======================================================================
0142 void FindVertex (const TheArc& A,
0143 const Handle(TheTopolTool)& Domain,
0144 TheFunction& Func,
0145 IntStart_SequenceOfPathPoint& pnt,
0146 const Standard_Real Toler)
0147 {
0148
0149 // Find the vertex of the arc A restriction solutions. It stores
0150 // Vertex in the list solutions pnt.
0151
0152
0153 TheVertex vtx;
0154 Standard_Real param,valf;
0155 Standard_Integer itemp;
0156
0157 Domain->Initialize(A);
0158 Domain->InitVertexIterator();
0159 while (Domain->MoreVertex()) {
0160 vtx = Domain->Vertex();
0161 param = TheSOBTool::Parameter(vtx,A);
0162
0163 // Evaluate the function and look compared to tolerance of the
0164 // Vertex. If distance <= tolerance then add a vertex to the list of solutions.
0165 // The arc is already assumed in the load function.
0166
0167 Func.Value(param,valf);
0168 if (Abs(valf) <= Toler) {
0169 itemp = Func.GetStateNumber();
0170 pnt.Append(IntStart_ThePathPoint(Func.Valpoint(itemp),Toler, vtx,A,param));
0171 // Solution is added
0172 }
0173 Domain->NextVertex();
0174 }
0175 }
0176
0177 Standard_Boolean IsDegenerated(const Handle(Adaptor3d_CurveOnSurface)& theCurve)
0178 {
0179 if (theCurve->GetType() == GeomAbs_Circle)
0180 {
0181 gp_Circ aCirc = theCurve->Circle();
0182 if (aCirc.Radius() <= Precision::Confusion())
0183 return Standard_True;
0184 }
0185 return Standard_False;
0186 }
0187
0188 Standard_Boolean IsDegenerated(const IntSurf_Quadric& theQuadric)
0189 {
0190 GeomAbs_SurfaceType TypeQuad = theQuadric.TypeQuadric();
0191 if (TypeQuad == GeomAbs_Cone)
0192 {
0193 gp_Cone aCone = theQuadric.Cone();
0194 Standard_Real aSemiAngle = Abs(aCone.SemiAngle());
0195 if (aSemiAngle < 0.02 || aSemiAngle > 1.55)
0196 return Standard_True;
0197 }
0198 return Standard_False;
0199 }
0200
0201 class SolInfo
0202 {
0203 public:
0204 SolInfo() : myMathIndex(-1), myValue(RealLast())
0205 {
0206 }
0207
0208 void Init(const math_FunctionAllRoots& theSolution, const Standard_Integer theIndex)
0209 {
0210 myMathIndex = theIndex;
0211 myValue = theSolution.GetPoint(theIndex);
0212 }
0213
0214 void Init(const IntCurveSurface_HInter& theSolution, const Standard_Integer theIndex)
0215 {
0216 myMathIndex = theIndex;
0217 myValue = theSolution.Point(theIndex).W();
0218 }
0219
0220 Standard_Real Value() const
0221 {
0222 return myValue;
0223 }
0224
0225 Standard_Integer Index() const
0226 {
0227 return myMathIndex;
0228 }
0229
0230 bool operator>(const SolInfo& theOther) const
0231 {
0232 return myValue > theOther.myValue;
0233 }
0234
0235 bool operator<(const SolInfo& theOther) const
0236 {
0237 return myValue < theOther.myValue;
0238 }
0239
0240 bool operator==(const SolInfo& theOther) const
0241 {
0242 return myValue == theOther.myValue;
0243 }
0244
0245 Standard_Real& ChangeValue()
0246 {
0247 return myValue;
0248 }
0249
0250 private:
0251 Standard_Integer myMathIndex;
0252 Standard_Real myValue;
0253 };
0254
0255 static
0256 void BoundedArc (const TheArc& A,
0257 const Handle(TheTopolTool)& Domain,
0258 const Standard_Real Pdeb,
0259 const Standard_Real Pfin,
0260 TheFunction& Func,
0261 IntStart_SequenceOfPathPoint& pnt,
0262 IntStart_SequenceOfSegment& seg,
0263 const Standard_Real TolBoundary,
0264 const Standard_Real TolTangency,
0265 Standard_Boolean& Arcsol,
0266 const Standard_Boolean RecheckOnRegularity)
0267 {
0268 // Recherche des points solutions et des bouts d arc solution sur un arc donne.
0269 // On utilise la fonction math_FunctionAllRoots. Ne convient donc que pour
0270 // des arcs ayant un point debut et un point de fin (intervalle ferme de
0271 // parametrage).
0272
0273 Standard_Integer i, Nbi = 0, Nbp = 0;
0274
0275 gp_Pnt ptdeb,ptfin;
0276 Standard_Real pardeb = 0., parfin = 0.;
0277 Standard_Integer ideb,ifin,range,ranged,rangef;
0278
0279 // Creer l echantillonage (math_FunctionSample ou classe heritant)
0280 // Appel a math_FunctionAllRoots
0281
0282 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0283 //@@@ La Tolerance est asociee a l arc ( Incoherence avec le cheminement )
0284 //@@@ ( EpsX ~ 1e-5 et ResolutionU et V ~ 1e-9 )
0285 //@@@ le vertex trouve ici n'est pas retrouve comme point d arret d une
0286 //@@@ ligne de cheminement
0287 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0288 Standard_Real EpsX = 1.e-10;
0289 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0290 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0291 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0292
0293 // Standard_Integer NbEchant = TheSOBTool::NbSamplesOnArc(A);
0294 Standard_Integer NbEchant = Func.NbSamples();
0295 if(NbEchant<100) NbEchant = 100; //-- lbr le 22 Avril 96
0296 //-- Toujours des pbs
0297
0298 //-- Modif 24 Aout 93 -----------------------------
0299 Standard_Real nTolTangency = TolTangency;
0300 if((Pfin - Pdeb) < (TolTangency*10.0)) {
0301 nTolTangency=(Pfin-Pdeb)*0.1;
0302 }
0303 if(EpsX>(nTolTangency+nTolTangency)) {
0304 EpsX = nTolTangency * 0.1;
0305 }
0306
0307 //--------------------------------------------------
0308 //-- Plante avec un edge avec 2 Samples
0309 //-- dont les extremites son solutions (f=0)
0310 //-- et ou la derivee est nulle
0311 //-- Exemple : un segment diametre d une sphere
0312 //-- if(NbEchant<3) NbEchant = 3; //-- lbr le 19 Avril 95
0313 //--------------------------------------------------
0314 Standard_Real para=0,dist,maxdist;
0315
0316 //-------------------------------------------------------------- REJECTIONS le 15 oct 98
0317 Standard_Boolean Rejection=Standard_True;
0318 Standard_Real maxdr,maxr,minr,ur,dur;
0319 minr=RealLast();
0320 maxr=-minr;
0321 maxdr=-minr;
0322 dur=(Pfin-Pdeb)*0.2;
0323 for(i=1,ur=Pdeb;i<=6;i++) {
0324 Standard_Real F,D;
0325 if(Func.Values(ur,F,D)) {
0326 Standard_Real lminr,lmaxr;
0327 if(D<0.0) D=-D;
0328 D*=dur+dur;
0329 if(D>maxdr) maxdr=D;
0330 lminr=F-D;
0331 lmaxr=F+D;
0332 if(lminr<minr) minr=lminr;
0333 if(lmaxr>maxr) maxr=lmaxr;
0334 if(minr<0.0 && maxr>0.0) {
0335 Rejection=Standard_False;
0336 break;
0337 }
0338 }
0339 ur+=dur;
0340 }
0341 if(Rejection)
0342 {
0343 dur=0.001+maxdr+(maxr-minr)*0.1;
0344 minr-=dur;
0345 maxr+=dur;
0346 if(minr<0.0 && maxr>0.0) {
0347 Rejection=Standard_False;
0348 }
0349 }
0350
0351 Arcsol=Standard_False;
0352
0353 if(Rejection==Standard_False)
0354 {
0355 const IntSurf_Quadric& aQuadric = Func.Quadric();
0356 GeomAbs_SurfaceType TypeQuad = aQuadric.TypeQuadric();
0357 GeomAbs_CurveType TypeConS = GeomAbs_OtherCurve;
0358
0359 IntCurveSurface_HInter IntCS;
0360 Standard_Boolean IsIntCSdone = Standard_False;
0361 TColStd_SequenceOfReal Params;
0362
0363 std::unique_ptr<math_FunctionAllRoots> pSol;
0364
0365 math_FunctionSample Echant(Pdeb,Pfin,NbEchant);
0366
0367 Standard_Boolean aelargir=Standard_True;
0368 //modified by NIZNHY-PKV Thu Apr 12 09:25:19 2001 f
0369 //
0370 //maxdist = 100.0*TolBoundary;
0371 maxdist = TolBoundary+TolTangency;
0372 //
0373 //modified by NIZNHY-PKV Thu Apr 12 09:25:23 2001 t
0374 for(i=1; i<=NbEchant && aelargir;i++) {
0375 Standard_Real u = Echant.GetParameter(i);
0376 if(Func.Value(u,dist)) {
0377 if(dist>maxdist || -dist>maxdist) {
0378 aelargir=Standard_False;
0379 }
0380 }
0381 }
0382 if(!(aelargir && maxdist<0.01)) {
0383 maxdist = TolBoundary;
0384 }
0385
0386 if (TypeQuad != GeomAbs_OtherSurface) //intersection of boundary curve and quadric surface
0387 {
0388 //Exact solution
0389 Handle(Adaptor3d_Surface) aSurf = Func.Surface();
0390 Adaptor3d_CurveOnSurface ConS(A, aSurf);
0391 TypeConS = ConS.GetType();
0392 #ifdef OCCT_DEBUG
0393 Handle(Geom_Curve) CurveConS;
0394 switch(TypeConS)
0395 {
0396 case GeomAbs_Line:
0397 {
0398 CurveConS = new Geom_Line(ConS.Line());
0399 break;
0400 }
0401 case GeomAbs_Circle:
0402 {
0403 CurveConS = new Geom_Circle(ConS.Circle());
0404 break;
0405 }
0406 case GeomAbs_Ellipse:
0407 {
0408 CurveConS = new Geom_Ellipse(ConS.Ellipse());
0409 break;
0410 }
0411 case GeomAbs_Hyperbola:
0412 {
0413 CurveConS = new Geom_Hyperbola(ConS.Hyperbola());
0414 break;
0415 }
0416 case GeomAbs_Parabola:
0417 {
0418 CurveConS = new Geom_Parabola(ConS.Parabola());
0419 break;
0420 }
0421 case GeomAbs_BezierCurve:
0422 {
0423 CurveConS = ConS.Bezier();
0424 break;
0425 }
0426 case GeomAbs_BSplineCurve:
0427 {
0428 CurveConS = ConS.BSpline();
0429 break;
0430 }
0431 default:
0432 {
0433 Standard_Real MaxDeviation, AverageDeviation;
0434 GeomLib::BuildCurve3d(1.e-5, ConS, ConS.FirstParameter(), ConS.LastParameter(),
0435 CurveConS, MaxDeviation, AverageDeviation);
0436 break;
0437 }
0438 }
0439 #endif
0440 Handle(Adaptor3d_CurveOnSurface) HConS = new Adaptor3d_CurveOnSurface(ConS);
0441 Handle(Geom_Surface) QuadSurf;
0442 switch (TypeQuad)
0443 {
0444 case GeomAbs_Plane:
0445 {
0446 QuadSurf = new Geom_Plane(aQuadric.Plane());
0447 break;
0448 }
0449 case GeomAbs_Cylinder:
0450 {
0451 QuadSurf = new Geom_CylindricalSurface(aQuadric.Cylinder());
0452 break;
0453 }
0454 case GeomAbs_Cone:
0455 {
0456 QuadSurf = new Geom_ConicalSurface(aQuadric.Cone());
0457 break;
0458 }
0459 case GeomAbs_Sphere:
0460 {
0461 QuadSurf = new Geom_SphericalSurface(aQuadric.Sphere());
0462 break;
0463 }
0464 case GeomAbs_Torus:
0465 {
0466 QuadSurf = new Geom_ToroidalSurface(aQuadric.Torus());
0467 break;
0468 }
0469 default:
0470 break;
0471 }
0472 Handle(GeomAdaptor_Surface) GAHsurf = new GeomAdaptor_Surface(QuadSurf);
0473
0474 if ((TypeConS == GeomAbs_Line ||
0475 TypeConS == GeomAbs_Circle ||
0476 TypeConS == GeomAbs_Ellipse ||
0477 TypeConS == GeomAbs_Parabola ||
0478 TypeConS == GeomAbs_Hyperbola) &&
0479 TypeQuad != GeomAbs_Torus &&
0480 !IsDegenerated(HConS) &&
0481 !IsDegenerated(aQuadric))
0482 {
0483 //exact intersection for only canonic curves and real quadric surfaces
0484 IntCS.Perform(HConS, GAHsurf);
0485 }
0486
0487 IsIntCSdone = IntCS.IsDone();
0488 if (IsIntCSdone)
0489 {
0490 Nbp = IntCS.NbPoints();
0491 Nbi = IntCS.NbSegments();
0492 }
0493 //If we have not got intersection, it may be touch with some tolerance,
0494 //need to be checked
0495 if (Nbp == 0 && Nbi == 0)
0496 IsIntCSdone = Standard_False;
0497
0498 } //if (TypeQuad != GeomAbs_OtherSurface) - intersection of boundary curve and quadric surface
0499
0500 if (!IsIntCSdone)
0501 {
0502 pSol.reset(new math_FunctionAllRoots(Func,Echant,EpsX,maxdist,maxdist)); //-- TolBoundary,nTolTangency);
0503
0504 if (!pSol->IsDone()) {throw Standard_Failure();}
0505
0506 Nbp=pSol->NbPoints();
0507 }
0508 //
0509 //jgv: build solution on the whole boundary
0510 if (RecheckOnRegularity && Nbp > 0 && IsRegularity(A, Domain))
0511 {
0512 //Standard_Real theTol = Domain->MaxTolerance(A);
0513 //theTol += theTol;
0514 Standard_Real theTol = 5.e-4;
0515 math_FunctionAllRoots SolAgain(Func,Echant,EpsX,theTol,theTol); //-- TolBoundary,nTolTangency);
0516
0517 if (!SolAgain.IsDone()) {throw Standard_Failure();}
0518
0519 Standard_Integer Nbi_again = SolAgain.NbIntervals();
0520
0521 if (Nbi_again > 0)
0522 {
0523 Standard_Integer NbSamples = 10;
0524 Standard_Real delta = (Pfin - Pdeb)/NbSamples;
0525 Standard_Real GlobalTol = theTol*10;
0526 Standard_Boolean SolOnBoundary = Standard_True;
0527 for (i = 0; i <= NbSamples; i++)
0528 {
0529 Standard_Real aParam = Pdeb + i*delta;
0530 Standard_Real aValue;
0531 Func.Value(aParam, aValue);
0532 if (Abs(aValue) > GlobalTol)
0533 {
0534 SolOnBoundary = Standard_False;
0535 break;
0536 }
0537 }
0538
0539 if (SolOnBoundary)
0540 {
0541 for (i = 1; i <= Nbi_again; i++)
0542 {
0543 IntStart_TheSegment newseg;
0544 newseg.SetValue(A);
0545 // Recuperer point debut et fin, et leur parametre.
0546 SolAgain.GetInterval(i,pardeb,parfin);
0547
0548 if (Abs(pardeb - Pdeb) <= Precision::PConfusion())
0549 pardeb = Pdeb;
0550 if (Abs(parfin - Pfin) <= Precision::PConfusion())
0551 parfin = Pfin;
0552
0553 SolAgain.GetIntervalState(i,ideb,ifin);
0554
0555 //-- cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx : i= "<<i<<" ParDeb:"<<pardeb<<" ParFin:"<<parfin<<endl;
0556
0557 ptdeb=Func.Valpoint(ideb);
0558 ptfin=Func.Valpoint(ifin);
0559
0560 PointProcess(ptdeb,pardeb,A,Domain,pnt,theTol,ranged);
0561 newseg.SetLimitPoint(pnt.Value(ranged),Standard_True);
0562 PointProcess(ptfin,parfin,A,Domain,pnt,theTol,rangef);
0563 newseg.SetLimitPoint(pnt.Value(rangef),Standard_False);
0564 seg.Append(newseg);
0565 }
0566 Arcsol=Standard_True;
0567 return;
0568 }
0569 }
0570 } //if (RecheckOnRegularity && Nbp > 0 && IsRegularity(A, Domain))
0571 ////////////////////////////////////////////
0572
0573 //-- detection du cas ou la fonction est quasi tangente et que les
0574 //-- zeros sont quasi confondus.
0575 //-- Dans ce cas on prend le point "milieu"
0576 //-- On suppose que les solutions sont triees.
0577
0578 if(Nbp) {
0579 NCollection_Array1<SolInfo> aSI(1, Nbp);
0580
0581 for(i=1;i<=Nbp;i++)
0582 {
0583 if (IsIntCSdone)
0584 aSI(i).Init(IntCS, i);
0585 else
0586 aSI(i).Init(*pSol, i);
0587 }
0588
0589 std::sort(aSI.begin(), aSI.end());
0590
0591 //modified by NIZNHY-PKV Wed Mar 21 18:34:18 2001 f
0592 //////////////////////////////////////////////////////////
0593 // The treatment of the situation when line(arc) that is
0594 // tangent to cylinder(domain).
0595 // We should have only one solution i.e Nbp=1. Ok?
0596 // But we have 2,3,.. solutions. That is wrong ersult.
0597 // The TreatLC(...) function is dedicated to solve the pb.
0598 // PKV Fri Mar 23 12:17:29 2001
0599
0600 Standard_Integer ip = TreatLC (A, Domain, aQuadric, TolBoundary, pnt);
0601 if (ip) {
0602 //////////////////////////////////////////////////////////
0603 //modified by NIZNHY-PKV Wed Mar 21 18:34:23 2001 t
0604 //
0605 // Using of old usual way proposed by Laurent
0606 //
0607 for(i=1;i<Nbp;i++) {
0608 Standard_Real parap1 = aSI(i + 1).Value();
0609 para = aSI(i).Value();
0610
0611 Standard_Real param=(para+parap1)*0.5;
0612 Standard_Real yf = 0.0;
0613 Standard_Real ym = 0.0;
0614 Standard_Real yl = 0.0;
0615 if(Func.Value(param,ym) && Abs(ym) < maxdist) {
0616 Standard_Real sm = Sign(1., ym);
0617 Standard_Boolean aTang = Func.Value(para,yf) && Func.Value(parap1,yl);
0618 if (aTang) {
0619 //Line can be tangent surface if all distances less then maxdist
0620 aTang = aTang && Abs(yf) < maxdist && Abs(yl) < maxdist;
0621 }
0622 if (aTang && IsIntCSdone && TypeConS == GeomAbs_Line) {
0623 //Interval is got by exact intersection
0624 //Line can be tangent if all points are on the same side of surface
0625 //it means that signs of all distances are the same
0626 Standard_Real sf = Sign(1., yf), sl = Sign(1., yl);
0627 aTang = aTang && (sm == sf) && (sm == sl);
0628 }
0629 if(aTang) {
0630 // Modified by skv - Tue Aug 31 12:13:51 2004 OCC569 Begin
0631 // Consider this interval as tangent one. Treat it to find
0632 // parameter with the lowest function value.
0633 // Compute the number of nodes.
0634 Standard_Real aTol = TolBoundary*1000.0;
0635 if(aTol > 0.001)
0636 aTol = 0.001;
0637
0638 // fix floating point exception 569, chl-922-e9
0639 parap1 = (Abs(parap1) < 1.e9) ? parap1 : ((parap1 >= 0.) ? 1.e9 : -1.e9);
0640 para = (Abs(para) < 1.e9) ? para : ((para >= 0.) ? 1.e9 : -1.e9);
0641
0642 Standard_Integer aNbNodes = RealToInt(Ceiling((parap1 - para)/aTol));
0643
0644 Standard_Real aVal = RealLast();
0645 Standard_Real aValMax = 0.;
0646 //Standard_Integer aNbNodes = 23;
0647 Standard_Real aDelta = (parap1 - para)/(aNbNodes + 1.);
0648 Standard_Integer ii;
0649 Standard_Real aCurPar;
0650 Standard_Real aCurVal;
0651
0652 for (ii = 0; ii <= aNbNodes + 1; ii++) {
0653 aCurPar = (ii < aNbNodes + 1) ? para + ii*aDelta : parap1;
0654
0655 if (Func.Value(aCurPar, aCurVal)) {
0656 Standard_Real anAbsVal = Abs(aCurVal);
0657 if (anAbsVal < aVal) {
0658 aVal = anAbsVal;
0659 param = aCurPar;
0660 }
0661 if (anAbsVal > aValMax)
0662 {
0663 aValMax = anAbsVal;
0664 }
0665 }
0666 }
0667 // At last, interval got by exact intersection can be considered as tangent if
0668 // minimal distance is inside interval and
0669 // minimal and maximal values are almost the same
0670 if (IsIntCSdone && aNbNodes > 1) {
0671 aTang = Abs(param - para) > EpsX && Abs(parap1 - param) > EpsX &&
0672 0.01*aValMax <= aVal;
0673 }
0674 if (aTang)
0675 {
0676 aSI(i).ChangeValue() = Pdeb - 1;
0677 aSI(i + 1).ChangeValue() = param;
0678 }
0679 }
0680 }
0681 }
0682
0683 for (i=1; i<=Nbp; i++) {
0684 para = aSI(i).Value();
0685 if((para-Pdeb)<EpsX || (Pfin-para)<EpsX)
0686 continue;
0687
0688 if(!Func.Value(para,dist))
0689 continue;
0690
0691 dist = Abs(dist);
0692
0693 Standard_Integer anIndx = -1;
0694 //const Standard_Real aParam = Sol->GetPoint(aSI(i).Index());
0695 const Standard_Real aParam = aSI(i).Value();
0696 if (dist < maxdist)
0697 {
0698 if (!IsIntCSdone &&
0699 (Abs(aParam - Pdeb) <= Precision::PConfusion() || Abs(aParam - Pfin) <= Precision::PConfusion()))
0700 {
0701 anIndx = pSol->GetPointState(aSI(i).Index());
0702 }
0703 }
0704
0705 gp_Pnt aPnt(anIndx < 0 ? Func.LastComputedPoint() : Func.Valpoint(anIndx));
0706
0707 if (dist > 0.1*Precision::Confusion())
0708 {
0709 //Precise found points. It results in following:
0710 // 1. Make the vertex nearer to the intersection line
0711 // (see description to issue #27252 in order to
0712 // understand necessity).
0713 // 2. Merge two near vertices to single point.
0714
0715 //All members in TabSol array has already been sorted in increase order.
0716 //Now, we limit precise boundaries in order to avoid changing this order.
0717 const Standard_Real aFPar = (i == 1) ? Pdeb : (para + aSI(i - 1).Value()) / 2.0;
0718 const Standard_Real aLPar = (i == Nbp) ? Pfin : (para + aSI(i + 1).Value()) / 2.0;
0719
0720 MinFunction aNewFunc(Func);
0721 math_BrentMinimum aMin(Precision::Confusion());
0722
0723 aMin.Perform(aNewFunc, aFPar, para, aLPar);
0724 if(aMin.IsDone())
0725 {
0726 para = aMin.Location();
0727 const gp_Pnt2d aP2d(A->Value(para));
0728 aPnt = Func.Surface()->Value(aP2d.X(), aP2d.Y());
0729 }
0730 }
0731
0732 PointProcess(aPnt, para, A, Domain, pnt, TolBoundary, range);
0733 }
0734 }// end of if(ip)
0735 } // end of if(Nbp)
0736
0737 // Pour chaque intervalle trouve faire
0738 // Traiter les extremites comme des points
0739 // Ajouter intervalle dans la liste des segments
0740
0741 if (!IsIntCSdone)
0742 Nbi = pSol->NbIntervals();
0743
0744 if (!RecheckOnRegularity && Nbp) {
0745 //--cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx :Nbp>0 0 <- Nbi "<<Nbi<<endl;
0746 Nbi=0;
0747 }
0748
0749 //-- cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx : Nbi : "<<Nbi<<endl;
0750
0751 for (i=1; i<=Nbi; i++) {
0752 IntStart_TheSegment newseg;
0753 newseg.SetValue(A);
0754 // Recuperer point debut et fin, et leur parametre.
0755 if (IsIntCSdone)
0756 {
0757 IntCurveSurface_IntersectionSegment IntSeg = IntCS.Segment(i);
0758 IntCurveSurface_IntersectionPoint End1 = IntSeg.FirstPoint();
0759 IntCurveSurface_IntersectionPoint End2 = IntSeg.SecondPoint();
0760 pardeb = End1.W();
0761 parfin = End2.W();
0762 ptdeb = End1.Pnt();
0763 ptfin = End2.Pnt();
0764 }
0765 else
0766 {
0767 pSol->GetInterval(i,pardeb,parfin);
0768 pSol->GetIntervalState(i,ideb,ifin);
0769
0770 //-- cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx : i= "<<i<<" ParDeb:"<<pardeb<<" ParFin:"<<parfin<<endl;
0771
0772 ptdeb=Func.Valpoint(ideb);
0773 ptfin=Func.Valpoint(ifin);
0774 }
0775
0776 PointProcess(ptdeb,pardeb,A,Domain,pnt,TolBoundary,ranged);
0777 newseg.SetLimitPoint(pnt.Value(ranged),Standard_True);
0778 PointProcess(ptfin,parfin,A,Domain,pnt,TolBoundary,rangef);
0779 newseg.SetLimitPoint(pnt.Value(rangef),Standard_False);
0780 seg.Append(newseg);
0781 }
0782
0783 if (Nbi==1) {
0784 if((Abs(pardeb - Pdeb) < Precision::PConfusion()) &&
0785 (Abs(parfin - Pfin) < Precision::PConfusion()))
0786 {
0787 Arcsol=Standard_True;
0788 }
0789 }
0790 }
0791 }
0792
0793 //=======================================================================
0794 //function : ComputeBoundsfromInfinite
0795 //purpose :
0796 //=======================================================================
0797 // - PROVISIONAL - TEMPORARY - NOT GOOD - NYI - TO DO
0798 // - Temporary - temporary - not good - nyi - to do
0799 void ComputeBoundsfromInfinite(TheFunction& Func,
0800 Standard_Real& PDeb,
0801 Standard_Real& PFin,
0802 Standard_Integer& NbEchant)
0803 {
0804
0805 // - We are looking for parameters for start and end of the arc (2d curve)
0806 // - Infinity, a way to intersect the quadric with a portion of arc
0807 // - Finished.
0808 //
0809 // - The quadric is a plane, a cylinder, a cone and a sphere.
0810 // - Idea: We take any point on the arc and the fact grow
0811 // - Terminals to the signed distance function values or is likely
0812 // - S cancel.
0813 //
0814 // - WARNING: The following calculations provide a very estimated coarse parameters.
0815 // - This avoids the raises and allows a case of Boxes
0816 // - Inifinies walk. It will take this code
0817 // - With curve surface intersections.
0818
0819 NbEchant = 100;
0820
0821 Standard_Real U0 = 0.0;
0822 Standard_Real dU = 0.001;
0823 Standard_Real Dist0,Dist1;
0824
0825 Func.Value(U0 , Dist0);
0826 Func.Value(U0+dU, Dist1);
0827 Standard_Real dDist = Dist1 - Dist0;
0828 if(dDist) {
0829 U0 -= dU*Dist0 / dDist;
0830 PDeb = PFin = U0;
0831 Standard_Real Umin = U0 - 1e5;
0832 Func.Value(Umin , Dist0);
0833 Func.Value(Umin+dU, Dist1);
0834 dDist = Dist1-Dist0;
0835 if(dDist) {
0836 Umin -= dU*Dist0 / dDist;
0837 }
0838 else {
0839 Umin-=10.0;
0840 }
0841 Standard_Real Umax = U0 + 1e8;
0842 Func.Value(Umax , Dist0);
0843 Func.Value(Umax+dU, Dist1);
0844 dDist = Dist1-Dist0;
0845 if(dDist) {
0846 Umax -= dU*Dist0 / dDist;
0847 }
0848 else {
0849 Umax+=10.0;
0850 }
0851 if(Umin>U0) { Umin=U0-10.0; }
0852 if(Umax<U0) { Umax=U0+10.0; }
0853
0854 PFin = Umax + 10. * (Umax - Umin);
0855 PDeb = Umin - 10. * (Umax - Umin);
0856 }
0857 else {
0858 //-- Possibilite de Arc totalement inclu ds Quad
0859 PDeb = 1e10;
0860 PFin = -1e10;
0861 }
0862 }
0863
0864 //=======================================================================
0865 //function : PointProcess
0866 //purpose :
0867 //=======================================================================
0868 void PointProcess (const gp_Pnt& Pt,
0869 const Standard_Real Para,
0870 const TheArc& A,
0871 const Handle(TheTopolTool)& Domain,
0872 IntStart_SequenceOfPathPoint& pnt,
0873 const Standard_Real Tol,
0874 Standard_Integer& Range)
0875 {
0876
0877 // Check to see if a solution point is coincident with a vertex.
0878 // If confused, you should find this vertex in the list of
0879 // Start. It then returns the position of this point in the list pnt.
0880 // Otherwise, add the point in the list.
0881
0882 Standard_Integer k;
0883 Standard_Boolean found,goon;
0884 Standard_Real dist,toler;
0885
0886 Standard_Integer Nbsol = pnt.Length();
0887 TheVertex vtx;
0888 IntStart_ThePathPoint ptsol;
0889
0890 Domain->Initialize(A);
0891 Domain->InitVertexIterator();
0892 found = Standard_False;
0893 goon = Domain->MoreVertex();
0894 while (goon) {
0895 vtx = Domain->Vertex();
0896 dist= Abs(Para-TheSOBTool::Parameter(vtx,A));
0897 toler = TheSOBTool::Tolerance(vtx,A);
0898 #ifdef OCCT_DEBUG
0899 if(toler>0.1) {
0900 std::cout<<"IntStart_SearchOnBoundaries_1.gxx : ** WARNING ** Tol Vertex="<<toler<<std::endl;
0901 std::cout<<" Ou Edge degenere Ou Kro pointu"<<std::endl;
0902 if(toler>10000) toler=1e-7;
0903 }
0904 #endif
0905
0906 if (dist <= toler) {
0907 // Locate the vertex in the list of solutions
0908 k=1;
0909 found = (k>Nbsol);
0910 while (!found) {
0911 ptsol = pnt.Value(k);
0912 if (!ptsol.IsNew()) {
0913 //jag 940608 if (ptsol.Vertex() == vtx && ptsol.Arc() == A) {
0914 if (Domain->Identical(ptsol.Vertex(),vtx) &&
0915 ptsol.Arc() == A &&
0916 Abs(ptsol.Parameter()-Para) <= toler) {
0917 found=Standard_True;
0918 }
0919 else {
0920 k=k+1;
0921 found=(k>Nbsol);
0922 }
0923 }
0924 else {
0925 k=k+1;
0926 found=(k>Nbsol);
0927 }
0928 }
0929 if (k<=Nbsol) { // We find the vertex
0930 Range = k;
0931 }
0932 else { // Otherwise
0933 ptsol.SetValue(Pt,Tol,vtx,A,Para);
0934 pnt.Append(ptsol);
0935 Range = pnt.Length();
0936 }
0937 found = Standard_True;
0938 goon = Standard_False;
0939 }
0940 else {
0941 Domain->NextVertex();
0942 goon = Domain->MoreVertex();
0943 }
0944 }
0945
0946 if (!found) { // No one is falling on a vertex
0947 //jgv: do not add segment's extremities if they already exist
0948 Standard_Boolean found_internal = Standard_False;
0949 for (k = 1; k <= pnt.Length(); k++)
0950 {
0951 ptsol = pnt.Value(k);
0952 if (ptsol.Arc() != A ||
0953 !ptsol.IsNew()) //vertex
0954 continue;
0955 if (Abs(ptsol.Parameter()-Para) <= Precision::PConfusion())
0956 {
0957 found_internal = Standard_True;
0958 Range = k;
0959 }
0960 }
0961 /////////////////////////////////////////////////////////////
0962
0963 if (!found_internal)
0964 {
0965 Standard_Real TOL=Tol;
0966 TOL*=1000.0;
0967 //if(TOL>0.001) TOL=0.001;
0968 if(TOL>0.005) TOL=0.005; //#24643
0969
0970 ptsol.SetValue(Pt,TOL,A,Para);
0971 pnt.Append(ptsol);
0972 Range = pnt.Length();
0973 }
0974 }
0975 }
0976
0977 //=======================================================================
0978 //function : IsRegularity
0979 //purpose :
0980 //=======================================================================
0981 Standard_Boolean IsRegularity(const TheArc& /*A*/,
0982 const Handle(TheTopolTool)& aDomain)
0983 {
0984 Standard_Address anEAddress=aDomain->Edge();
0985 if (anEAddress==NULL) {
0986 return Standard_False;
0987 }
0988
0989 TopoDS_Edge* anE=(TopoDS_Edge*)anEAddress;
0990
0991 return (BRep_Tool::HasContinuity(*anE));
0992 }
0993
0994 //=======================================================================
0995 //function : TreatLC
0996 //purpose :
0997 //=======================================================================
0998 Standard_Integer TreatLC (const TheArc& A,
0999 const Handle(TheTopolTool)& aDomain,
1000 const IntSurf_Quadric& aQuadric,
1001 const Standard_Real TolBoundary,
1002 IntStart_SequenceOfPathPoint& pnt)
1003 {
1004 Standard_Integer anExitCode=1, aNbExt;
1005
1006 Standard_Address anEAddress=aDomain->Edge();
1007 if (anEAddress==NULL) {
1008 return anExitCode;
1009 }
1010
1011 TopoDS_Edge* anE=(TopoDS_Edge*)anEAddress;
1012
1013 if (BRep_Tool::Degenerated(*anE)) {
1014 return anExitCode;
1015 }
1016
1017 GeomAbs_CurveType aTypeE;
1018 BRepAdaptor_Curve aBAC(*anE);
1019 aTypeE=aBAC.GetType();
1020
1021 if (aTypeE!=GeomAbs_Line) {
1022 return anExitCode;
1023 }
1024
1025 GeomAbs_SurfaceType aTypeS;
1026 aTypeS=aQuadric.TypeQuadric();
1027
1028 if (aTypeS!=GeomAbs_Cylinder) {
1029 return anExitCode;
1030 }
1031
1032 Standard_Real f, l, U1f, U1l, U2f, U2l, UEgde, TOL, aDist, aR, aRRel, Tol;
1033 Handle(Geom_Curve) aCEdge=BRep_Tool::Curve(*anE, f, l);
1034
1035 gp_Cylinder aCyl=aQuadric.Cylinder();
1036 const gp_Ax1& anAx1=aCyl.Axis();
1037 gp_Lin aLin(anAx1);
1038 Handle(Geom_Line) aCAxis=new Geom_Line (aLin);
1039 aR=aCyl.Radius();
1040
1041 U1f = aCAxis->FirstParameter();
1042 U1l = aCAxis->LastParameter();
1043
1044 U2f = aCEdge->FirstParameter();
1045 U2l = aCEdge->LastParameter();
1046
1047
1048 GeomAdaptor_Curve C1, C2;
1049
1050 C1.Load(aCAxis);
1051 C2.Load(aCEdge);
1052
1053 Tol = Precision::PConfusion();
1054
1055 Extrema_ExtCC anExtCC(C1, C2, U1f, U1l, U2f, U2l, Tol, Tol);
1056
1057 aNbExt=anExtCC.NbExt();
1058 if (aNbExt!=1) {
1059 return anExitCode;
1060 }
1061
1062 gp_Pnt P1,PEdge;
1063 Extrema_POnCurv PC1, PC2;
1064
1065 anExtCC.Points(1, PC1, PC2);
1066
1067 P1 =PC1.Value();
1068 PEdge=PC2.Value();
1069
1070 UEgde=PC2.Parameter();
1071
1072 aDist=PEdge.Distance(P1);
1073 aRRel=fabs(aDist-aR)/aR;
1074 if (aRRel > TolBoundary) {
1075 return anExitCode;
1076 }
1077
1078 if (UEgde < (f+TolBoundary) || UEgde > (l-TolBoundary)) {
1079 return anExitCode;
1080 }
1081 //
1082 // Do not wonder !
1083 // It was done as into PointProcess(...) function
1084 //printf("TreatLC()=> tangent line is found\n");
1085 TOL=1000.*TolBoundary;
1086 if(TOL>0.001) TOL=0.001;
1087
1088 IntStart_ThePathPoint ptsol;
1089 ptsol.SetValue(PEdge, TOL, A, UEgde);
1090 pnt.Append(ptsol);
1091
1092 anExitCode=0;
1093 return anExitCode;
1094
1095 }
1096
1097
1098 //=======================================================================
1099 //function : IntStart_SearchOnBoundaries::IntStart_SearchOnBoundaries
1100 //purpose :
1101 //=======================================================================
1102 IntStart_SearchOnBoundaries::IntStart_SearchOnBoundaries ()
1103 : done(Standard_False),
1104 all(Standard_False)
1105 {
1106 }
1107
1108 //=======================================================================
1109 //function : Perform
1110 //purpose :
1111 //=======================================================================
1112 void IntStart_SearchOnBoundaries::Perform (TheFunction& Func,
1113 const Handle(TheTopolTool)& Domain,
1114 const Standard_Real TolBoundary,
1115 const Standard_Real TolTangency,
1116 const Standard_Boolean RecheckOnRegularity)
1117 {
1118
1119 done = Standard_False;
1120 spnt.Clear();
1121 sseg.Clear();
1122
1123 Standard_Boolean Arcsol;
1124 Standard_Real PDeb,PFin, prm, tol;
1125 Standard_Integer i, nbknown, nbfound,index;
1126 gp_Pnt pt;
1127
1128 Domain->Init();
1129
1130 if (Domain->More()) {
1131 all = Standard_True;
1132 }
1133 else {
1134 all = Standard_False;
1135 }
1136
1137 while (Domain->More()) {
1138 TheArc A = Domain->Value();
1139 if (!TheSOBTool::HasBeenSeen(A)) {
1140 Func.Set(A);
1141 FindVertex(A,Domain,Func,spnt,TolBoundary);
1142 TheSOBTool::Bounds(A,PDeb,PFin);
1143 if(Precision::IsNegativeInfinite(PDeb) ||
1144 Precision::IsPositiveInfinite(PFin)) {
1145 Standard_Integer NbEchant;
1146 ComputeBoundsfromInfinite(Func,PDeb,PFin,NbEchant);
1147 }
1148 BoundedArc(A,Domain,PDeb,PFin,Func,spnt,sseg,TolBoundary,TolTangency,Arcsol,RecheckOnRegularity);
1149 all = (all && Arcsol);
1150 }
1151
1152 else {
1153 // as it seems we'll never be here, because
1154 // TheSOBTool::HasBeenSeen(A) always returns FALSE
1155 nbfound = spnt.Length();
1156
1157 // On recupere les points connus
1158 nbknown = TheSOBTool::NbPoints(A);
1159 for (i=1; i<=nbknown; i++) {
1160 TheSOBTool::Value(A,i,pt,tol,prm);
1161 if (TheSOBTool::IsVertex(A,i)) {
1162 TheVertex vtx;
1163 TheSOBTool::Vertex(A,i,vtx);
1164 spnt.Append(IntStart_ThePathPoint(pt,tol,vtx,A,prm));
1165 }
1166 else {
1167 spnt.Append(IntStart_ThePathPoint(pt,tol,A,prm));
1168 }
1169 }
1170 // On recupere les arcs solutions
1171 nbknown = TheSOBTool::NbSegments(A);
1172 for (i=1; i<=nbknown; i++) {
1173 IntStart_TheSegment newseg;
1174 newseg.SetValue(A);
1175 if (TheSOBTool::HasFirstPoint(A,i,index)) {
1176 newseg.SetLimitPoint(spnt.Value(nbfound+index),Standard_True);
1177 }
1178 if (TheSOBTool::HasLastPoint(A,i,index)) {
1179 newseg.SetLimitPoint(spnt.Value(nbfound+index),Standard_False);
1180 }
1181 sseg.Append(newseg);
1182 }
1183 all = (all& TheSOBTool::IsAllSolution(A));
1184 }
1185 Domain->Next();
1186 }
1187 done = Standard_True;
1188 }
