diff --git a/dox/user_guides/modeling_data/modeling_data.md b/dox/user_guides/modeling_data/modeling_data.md index 1e90d98938..6c153507aa 100644 --- a/dox/user_guides/modeling_data/modeling_data.md +++ b/dox/user_guides/modeling_data/modeling_data.md @@ -338,7 +338,8 @@ The GeomConvert package also provides the following: * a splitting algorithm, which determines the curves along which a BSpline surface should be cut in order to obtain patches with the same degree of continuity, * global functions to construct BSpline surfaces created by this splitting algorithm, or by other types of BSpline surface segmentation, * an algorithm, which converts a BSpline surface into a series of adjacent Bezier surfaces, - * an algorithm, which converts a grid of adjacent Bezier surfaces into a BSpline surface. + * an algorithm, which converts a grid of adjacent Bezier surfaces into a BSpline surface. + * algorithms that converts NURBS, Bezier and other general parametrized curves and surface into anaytical curves and surfaces. @subsection occt_modat_1_4 Points on Curves diff --git a/dox/user_guides/shape_healing/shape_healing.md b/dox/user_guides/shape_healing/shape_healing.md index 0f811679ed..fe479e25ec 100644 --- a/dox/user_guides/shape_healing/shape_healing.md +++ b/dox/user_guides/shape_healing/shape_healing.md @@ -1007,6 +1007,25 @@ Standard_Integer aNbOffsetSurfaces = aCheckContents.NbOffsetSurf(); Handle(TopTools_HSequenceOfShape) aSeqFaces = aCheckContents.OffsetSurfaceSec(); ~~~~ +@subsubsection occt_shg_3_2_4 Analysis of shape underlined geometry + +Class *ShapeAnalysis_CanonicalRecognition* provides tools that analyze geometry of shape and explore the possibility of converting geometry into a canonical form. +Canonical forms for curves are lines, circles and ellipses. +Canonical forms for surfaces are planar, cylindrical, conical and spherical surfaces. + +Recognition and converting into canonical form is performed according to maximal deviation criterium: maximal distance between initial and canonical geometrical objects must be less, than given value. + +Analysis of curves is allowed for following shapes: + * edge - algorithm checks 3d curve of edge + * wire - algorithm checks 3d curves of all edges in order to convert them in the same analytical curve + +Analysis of surfaces is allowed for following shapes: + * face - algorithm checks surface of face + * shell - algorithm checks surfaces of all faces in order to convert them in the same analytical surface + * edge - algorithm checks all surfaces that are shared by given edge in order convert one of them in analytical surface, which most close to the input sample surface. + * wire - the same as for edge, but algorithm checks all edges of wire in order to find analytical surface, which most close to the input sample surface. + + @section occt_shg_4 Upgrading Upgrading tools are intended for adaptation of shapes for better use by Open CASCADE Technology or for customization to particular needs, i.e. for export to another system. diff --git a/src/GeomConvert/FILES b/src/GeomConvert/FILES index 029d1cb1a2..06ae4fe8f8 100755 --- a/src/GeomConvert/FILES +++ b/src/GeomConvert/FILES @@ -20,3 +20,8 @@ GeomConvert_CompCurveToBSplineCurve.cxx GeomConvert_CompCurveToBSplineCurve.hxx GeomConvert_Units.cxx GeomConvert_Units.hxx +GeomConvert_CurveToAnaCurve.cxx +GeomConvert_CurveToAnaCurve.hxx +GeomConvert_SurfToAnaSurf.cxx +GeomConvert_SurfToAnaSurf.hxx +GeomConvert_ConvType.hxx diff --git a/src/GeomConvert/GeomConvert_ConvType.hxx b/src/GeomConvert/GeomConvert_ConvType.hxx new file mode 100644 index 0000000000..56592c7aeb --- /dev/null +++ b/src/GeomConvert/GeomConvert_ConvType.hxx @@ -0,0 +1,23 @@ +// Copyright (c) 2022 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _GeomConvert_ConvType_HeaderFile +#define _GeomConvert_ConvType_HeaderFile +enum GeomConvert_ConvType +{ + GeomConvert_Target, + GeomConvert_Simplest, + GeomConvert_MinGap +}; + +#endif // _GeomConvert_ConvType_HeaderFile diff --git a/src/GeomConvert/GeomConvert_CurveToAnaCurve.cxx b/src/GeomConvert/GeomConvert_CurveToAnaCurve.cxx new file mode 100644 index 0000000000..8e135e1edb --- /dev/null +++ b/src/GeomConvert/GeomConvert_CurveToAnaCurve.cxx @@ -0,0 +1,768 @@ +// Created: 2001-05-21 +// +// Copyright (c) 2001-2013 OPEN CASCADE SAS +// +// This file is part of commercial software by OPEN CASCADE SAS, +// furnished in accordance with the terms and conditions of the contract +// and with the inclusion of this copyright notice. +// This file or any part thereof may not be provided or otherwise +// made available to any third party. +// +// No ownership title to the software is transferred hereby. +// +// OPEN CASCADE SAS makes no representation or warranties with respect to the +// performance of this software, and specifically disclaims any responsibility +// for any damages, special or consequential, connected with its use. + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +GeomConvert_CurveToAnaCurve::GeomConvert_CurveToAnaCurve(): + myGap(Precision::Infinite()), + myConvType(GeomConvert_MinGap), + myTarget(GeomAbs_Line) +{ +} + +GeomConvert_CurveToAnaCurve::GeomConvert_CurveToAnaCurve(const Handle(Geom_Curve)& C) : + myGap(Precision::Infinite()), + myConvType(GeomConvert_MinGap), + myTarget(GeomAbs_Line) +{ + myCurve = C; +} + +void GeomConvert_CurveToAnaCurve::Init(const Handle(Geom_Curve)& C) +{ + myCurve = C; + myGap = Precision::Infinite(); +} + +//======================================================================= +//function : ConvertToAnalytical +//purpose : +//======================================================================= + +Standard_Boolean GeomConvert_CurveToAnaCurve::ConvertToAnalytical(const Standard_Real tol, + Handle(Geom_Curve)& theResultCurve, + const Standard_Real F, const Standard_Real L, + Standard_Real& NewF, Standard_Real& NewL) +{ + if(myCurve.IsNull()) + return Standard_False; + + Handle(Geom_Curve) aCurve = myCurve; + while (aCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) { + Handle(Geom_TrimmedCurve) aTrimmed = Handle(Geom_TrimmedCurve):: + DownCast(aCurve); + aCurve = aTrimmed->BasisCurve(); + } + + Handle(Geom_Curve) C = ComputeCurve(aCurve,tol,F, L, NewF, NewL, myGap, myConvType, myTarget); + + if(C.IsNull()) return Standard_False; + theResultCurve = C; + return Standard_True; +} + +//======================================================================= +//function : IsLinear +//purpose : +//======================================================================= + +Standard_Boolean GeomConvert_CurveToAnaCurve::IsLinear(const TColgp_Array1OfPnt& aPoles, + const Standard_Real tolerance, + Standard_Real& Deviation) +{ + Standard_Integer nbPoles = aPoles.Length(); + if(nbPoles < 2) + return Standard_False; + + Standard_Real dMax = 0; + Standard_Integer iMax1=0,iMax2=0; + + Standard_Integer i; + for(i = 1; i < nbPoles; i++) + for(Standard_Integer j = i+1; j <= nbPoles; j++) { + Standard_Real dist = aPoles(i).SquareDistance(aPoles(j)); + if(dist > dMax) { + dMax = dist; + iMax1 = i; + iMax2 = j; + } + } + + if (dMax < Precision::SquareConfusion()) + return Standard_False; + + Standard_Real tol2 = tolerance*tolerance; + gp_Vec avec (aPoles(iMax1),aPoles(iMax2)); gp_Dir adir (avec); gp_Lin alin (aPoles(iMax1),adir); + + Standard_Real aMax = 0.; + for(i = 1; i <= nbPoles; i++) { + Standard_Real dist = alin.SquareDistance(aPoles(i)); + if(dist > tol2) + return Standard_False; + if(dist > aMax) + aMax = dist; + } + Deviation = sqrt(aMax); + + return Standard_True; +} + +//======================================================================= +//function : GetLine +//purpose : +//======================================================================= + +gp_Lin GeomConvert_CurveToAnaCurve::GetLine(const gp_Pnt& P1, const gp_Pnt& P2, + Standard_Real& cf, Standard_Real& cl) +{ + gp_Vec avec(P1, P2); gp_Dir adir(avec); gp_Lin alin(P1, adir); + cf = ElCLib::Parameter(alin, P1); + cl = ElCLib::Parameter(alin, P2); + return alin; +} +//======================================================================= +//function : ComputeLine +//purpose : +//======================================================================= + +Handle(Geom_Line) GeomConvert_CurveToAnaCurve::ComputeLine (const Handle(Geom_Curve)& curve, + const Standard_Real tolerance, + const Standard_Real c1, const Standard_Real c2, + Standard_Real& cf, Standard_Real& cl, + Standard_Real& Deviation) +{ + Handle(Geom_Line) line; + if (curve.IsNull()) return line; + line = Handle(Geom_Line)::DownCast(curve); // qui sait + if (!line.IsNull()) { + cf = c1; + cl = c2; + Deviation = 0.; + return line; + } + + gp_Pnt P1 = curve->Value (c1); + gp_Pnt P2 = curve->Value (c2); + if(P1.SquareDistance(P2) < Precision::SquareConfusion()) + return line; + cf = c1; cl = c2; + + Handle(TColgp_HArray1OfPnt) Poles; + Standard_Integer nbPoles; + Handle(Geom_BSplineCurve) bsc = Handle(Geom_BSplineCurve)::DownCast(curve); + if (!bsc.IsNull()) { + nbPoles = bsc->NbPoles(); + Poles = new TColgp_HArray1OfPnt(1, nbPoles); + bsc->Poles(Poles->ChangeArray1()); + } + else + { + Handle(Geom_BezierCurve) bzc = Handle(Geom_BezierCurve)::DownCast(curve); + if (!bzc.IsNull()) { + nbPoles = bzc->NbPoles(); + Poles = new TColgp_HArray1OfPnt(1, nbPoles); + bzc->Poles(Poles->ChangeArray1()); + } + else + { + nbPoles = 23; + Poles = new TColgp_HArray1OfPnt(1, nbPoles); + Standard_Real dt = (c2 - c1) / (nbPoles - 1); + Poles->SetValue(1, P1); + Poles->SetValue(nbPoles, P2); + Standard_Integer i; + for (i = 2; i < nbPoles; ++i) + { + Poles->SetValue(i, curve->Value(c1 + (i - 1) * dt)); + } + } + } + if(!IsLinear(Poles->Array1(),tolerance,Deviation)) return line; // non + gp_Lin alin = GetLine (P1, P2, cf, cl); + line = new Geom_Line (alin); + return line; +} + +//======================================================================= +//function : GetCircle +//purpose : +//======================================================================= + +Standard_Boolean GeomConvert_CurveToAnaCurve::GetCircle (gp_Circ& crc, + const gp_Pnt& P0,const gp_Pnt& P1, const gp_Pnt& P2) +{ +// Control if points are not aligned (should be done by MakeCirc + Standard_Real aMaxCoord = Sqrt(Precision::Infinite()); + if (Abs(P0.X()) > aMaxCoord || Abs(P0.Y()) > aMaxCoord || Abs(P0.Z()) > aMaxCoord) + return Standard_False; + if (Abs(P1.X()) > aMaxCoord || Abs(P1.Y()) > aMaxCoord || Abs(P1.Z()) > aMaxCoord) + return Standard_False; + if (Abs(P2.X()) > aMaxCoord || Abs(P2.Y()) > aMaxCoord || Abs(P2.Z()) > aMaxCoord) + return Standard_False; + +// Building the circle + gce_MakeCirc mkc (P0,P1,P2); + if (!mkc.IsDone()) return Standard_False; + crc = mkc.Value(); + if (crc.Radius() < gp::Resolution()) return Standard_False; + // Recalage sur P0 + gp_Pnt PC = crc.Location(); + gp_Ax2 axe = crc.Position(); + gp_Vec VX (PC,P0); + axe.SetXDirection (VX); + crc.SetPosition (axe); + return Standard_True; +} + +//======================================================================= +//function : ComputeCircle +//purpose : +//======================================================================= + +Handle(Geom_Curve) GeomConvert_CurveToAnaCurve::ComputeCircle (const Handle(Geom_Curve)& c3d, + const Standard_Real tol, + const Standard_Real c1, const Standard_Real c2, + Standard_Real& cf, Standard_Real& cl, + Standard_Real& Deviation) +{ + if (c3d->IsKind (STANDARD_TYPE(Geom_Circle))) { + cf = c1; + cl = c2; + Deviation = 0.; + Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(c3d); + return aCirc; + } + + Handle(Geom_Circle) circ; + gp_Pnt P0,P1,P2; + Standard_Real ca = (c1+c1+c2) / 3; Standard_Real cb = (c1+c2+c2) / 3; + P0 = c3d->Value(c1); + P1 = c3d->Value(ca); + P2 = c3d->Value(cb); + + gp_Circ crc; + if (!GetCircle (crc,P0,P1,P2)) return circ; + + // Reste a controler que c est bien un cercle : prendre 20 points + Standard_Real du = (c2-c1)/20; + Standard_Integer i; + Standard_Real aMax = 0.; + for (i = 0; i <= 20; i ++) { + Standard_Real u = c1+(du*i); + gp_Pnt PP = c3d->Value(u); + Standard_Real dist = crc.Distance(PP); + if (dist > tol) return circ; // not done + if (dist > aMax) + aMax = dist; + } + Deviation = aMax; + + // defining the parameters + Standard_Real PI2 = 2 * M_PI; + + cf = ElCLib::Parameter (crc,c3d->Value (c1)); + cf = ElCLib::InPeriod(cf, 0., PI2); + + //first parameter should be closed to zero + + if(Abs(cf) < Precision::PConfusion() || Abs(PI2-cf) < Precision::PConfusion()) + cf = 0.; + + Standard_Real cm = ElCLib::Parameter (crc,c3d->Value ((c1+c2)/2.)); + cm = ElCLib::InPeriod(cm, cf, cf + PI2); + + cl = ElCLib::Parameter (crc,c3d->Value (c2)); + cl = ElCLib::InPeriod(cl, cm, cm + PI2); + + circ = new Geom_Circle (crc); + return circ; +} + +//======================================================================= +// Compute Ellipse +//======================================================================= + +//======================================================================= +//function : IsArrayPntPlanar +//purpose : +//======================================================================= + +static Standard_Boolean IsArrayPntPlanar(const Handle(TColgp_HArray1OfPnt)& HAP, + gp_Dir& Norm, const Standard_Real prec) +{ + Standard_Integer size = HAP->Length(); + if(size<3) + return Standard_False; + gp_Pnt P1 = HAP->Value(1); + gp_Pnt P2 = HAP->Value(2); + gp_Pnt P3 = HAP->Value(3); + Standard_Real dist1 = P1.Distance(P2); + Standard_Real dist2 = P1.Distance(P3); + if( dist13) { + for(i=4; i<=size; i++) { + gp_Pnt PN = HAP->Value(i); + dist1 = P1.Distance(PN); + if (dist1 < prec || Precision::IsInfinite(dist1)) + { + return Standard_False; + } + gp_Vec VN(P1,PN); + if(!NV.IsNormal(VN,prec)) + return Standard_False; + } + } + Norm = NV; + return Standard_True; +} + +//======================================================================= +//function : ConicdDefinition +//purpose : +//======================================================================= + +static Standard_Boolean ConicDefinition + ( const Standard_Real a, const Standard_Real b1, const Standard_Real c, + const Standard_Real d1, const Standard_Real e1, const Standard_Real f, + const Standard_Boolean IsParab, const Standard_Boolean IsEllip, + gp_Pnt& Center, gp_Dir& MainAxis, Standard_Real& Rmin, Standard_Real& Rmax ) +{ + Standard_Real Xcen = 0.,Ycen = 0., Xax = 0.,Yax = 0.; + Standard_Real b,d,e; + // conic : a*x2 + 2*b*x*y + c*y2 + 2*d*x + 2*e*y + f = 0. + //Equation (a,b,c,d,e,f); + b = b1/2.; d = d1/2.; e = e1/2.; // chgt de variable + + Standard_Real eps = 1.E-08; // ?? comme ComputedForm + + if (IsParab) { + + } + else { + // -> Conique a centre, cas general + // On utilise les Determinants des matrices : + // | a b d | + // gdet (3x3) = | b c e | et pdet (2X2) = | a b | + // | d e f | | b c | + + Standard_Real gdet = a*c*f + 2*b*d*e - c*d*d - a*e*e - b*b*f; + Standard_Real pdet = a*c - b*b; + + Xcen = (b*e - c*d) / pdet; + Ycen = (b*d - a*e) / pdet; + + Standard_Real term1 = a-c; + Standard_Real term2 = 2*b; + Standard_Real cos2t; + Standard_Real auxil; + + if (Abs(term2) <= eps && Abs(term1) <= eps) { + cos2t = 1.; + auxil = 0.; + } + else { + if (Abs(term1) < eps) + { + return Standard_False; + } + Standard_Real t2d = term2/term1; //skl 21.11.2001 + cos2t = 1./sqrt(1+t2d*t2d); + auxil = sqrt (term1*term1 + term2*term2); + } + + Standard_Real cost = sqrt ( (1+cos2t)/2. ); + Standard_Real sint = sqrt ( (1-cos2t)/2. ); + + Standard_Real aprim = (a+c+auxil)/2.; + Standard_Real cprim = (a+c-auxil)/2.; + + if (Abs(aprim) < gp::Resolution() || Abs(cprim) < gp::Resolution()) + return Standard_False; + + term1 = -gdet/(aprim*pdet); + term2 = -gdet/(cprim*pdet); + + if (IsEllip) { + Xax = cost; + Yax = sint; + Rmin = sqrt ( term1); + Rmax = sqrt ( term2); + if(RmaxIsKind (STANDARD_TYPE(Geom_Ellipse))) { + cf = c1; + cl = c2; + Deviation = 0.; + Handle(Geom_Ellipse) anElips = Handle(Geom_Ellipse)::DownCast(c3d); + return anElips; + } + + Handle(Geom_Curve) res; + Standard_Real prec = Precision::PConfusion(); + + Standard_Real AF,BF,CF,DF,EF,Q1,Q2,Q3,c2n; + Standard_Integer i; + + gp_Pnt PStart = c3d->Value(c1); + gp_Pnt PEnd = c3d->Value(c2); + + const Standard_Boolean IsClos = PStart.Distance(PEnd) < prec; + if (IsClos) + { + c2n=c2-(c2-c1)/5; + } + else + c2n=c2; + // + gp_XYZ aBC; + Handle(TColgp_HArray1OfPnt) AP = new TColgp_HArray1OfPnt(1,5); + AP->SetValue(1,PStart); + aBC += PStart.XYZ(); + Standard_Real dc=(c2n-c1)/4; + for (i = 1; i < 5; i++) + { + gp_Pnt aP = c3d->Value(c1 + dc*i); + AP->SetValue(i + 1, aP); + aBC += aP.XYZ(); + } + aBC /= 5; + aBC *= -1; + gp_Vec aTrans(aBC); + for (i = 1; i <= 5; ++i) + { + AP->ChangeValue(i).Translate(aTrans); + } + gp_Dir ndir; + if(!IsArrayPntPlanar(AP,ndir,prec)) + return res; + + if (Abs(ndir.X()) < gp::Resolution() && Abs(ndir.Y()) < gp::Resolution() + && Abs(ndir.Z()) < gp::Resolution()) + return res; + + gp_Ax3 AX(gp_Pnt(0,0,0),ndir); + gp_Trsf Tr; + Tr.SetTransformation(AX); + gp_Trsf Tr2 = Tr.Inverted(); + + math_Matrix Dt(1, 5, 1, 5); + math_Vector F(1, 5), Sl(1, 5); + + Standard_Real XN,YN,ZN = 0.; + gp_Pnt PT,PP; + for(i=1; i<=5; i++) { + PT = AP->Value(i).Transformed(Tr); + PT.Coord(XN,YN,ZN); + Dt(i, 1) = XN*XN; + Dt(i, 2) = XN*YN; + Dt(i, 3) = YN*YN; + Dt(i, 4) = XN; + Dt(i, 5) = YN; + F(i) = -1.; + } + + math_Gauss aSolver(Dt); + if (!aSolver.IsDone()) + return res; + + aSolver.Solve(F, Sl); + + AF=Sl(1); + BF=Sl(2); + CF=Sl(3); + DF=Sl(4); + EF=Sl(5); + + Q1=AF*CF+BF*EF*DF/4-CF*DF*DF/4-BF*BF/4-AF*EF*EF/4; + Q2=AF*CF-BF*BF/4; + Q3=AF+CF; + + Standard_Real Rmax, Rmin; + gp_Pnt Center; + gp_Dir MainAxis; + Standard_Boolean IsParab = Standard_False, IsEllip = Standard_False; + +if (Q2 > 0 && Q1*Q3 < 0) { + // ellipse + IsEllip = Standard_True; + if (ConicDefinition(AF, BF, CF, DF, EF, 1., IsParab, IsEllip, + Center, MainAxis, Rmin, Rmax)) { + // create ellipse + if (Rmax - Rmin < Precision::Confusion()) + { + return res; //really it is circle, which must be recognized in other method + } + aTrans *= -1; + Center.SetZ(ZN); + gp_Pnt NewCenter = Center.Transformed(Tr2); + gp_Pnt Ptmp(Center.X() + MainAxis.X() * 10, + Center.Y() + MainAxis.Y() * 10, + Center.Z() + MainAxis.Z() * 10); + gp_Pnt NewPtmp = Ptmp.Transformed(Tr2); + gp_Dir NewMainAxis(NewPtmp.X() - NewCenter.X(), + NewPtmp.Y() - NewCenter.Y(), + NewPtmp.Z() - NewCenter.Z()); + gp_Ax2 ax2(NewCenter, ndir, NewMainAxis); + + gp_Elips anEllipse(ax2, Rmax, Rmin); + anEllipse.Translate(aTrans); + Handle(Geom_Ellipse) gell = new Geom_Ellipse(anEllipse); + + // test for 20 points + Standard_Real param2 = 0; + dc = (c2 - c1) / 20; + for (i = 1; i <= 20; i++) { + PP = c3d->Value(c1 + i*dc); + Standard_Real aPar = ElCLib::Parameter(anEllipse, PP); + Standard_Real dist = gell->Value(aPar).Distance(PP); + if (dist > tol) return res; // not done + if (dist > param2) + param2 = dist; + } + + + Deviation = param2; + + Standard_Real PI2 = 2 * M_PI; + cf = ElCLib::Parameter(anEllipse, c3d->Value(c1)); + cf = ElCLib::InPeriod(cf, 0., PI2); + + //first parameter should be closed to zero + + if (Abs(cf) < Precision::PConfusion() || Abs(PI2 - cf) < Precision::PConfusion()) + cf = 0.; + + Standard_Real cm = ElCLib::Parameter(anEllipse, c3d->Value((c1 + c2) / 2.)); + cm = ElCLib::InPeriod(cm, cf, cf + PI2); + + cl = ElCLib::Parameter(anEllipse, c3d->Value(c2)); + cl = ElCLib::InPeriod(cl, cm, cm + PI2); + + res = gell; + } +} +/* +if (Q2 < 0 && Q1 != 0) { + // hyberbola +} + +if (Q2 == 0 && Q1 != 0) { + // parabola +} +*/ +return res; +} + + +//======================================================================= +//function : ComputeCurve +//purpose : +//======================================================================= + +Handle(Geom_Curve) GeomConvert_CurveToAnaCurve::ComputeCurve(const Handle(Geom_Curve)& theC3d, + const Standard_Real tolerance, + const Standard_Real c1, const Standard_Real c2, + Standard_Real& cf, Standard_Real& cl, + Standard_Real& theGap, + const GeomConvert_ConvType theConvType, const GeomAbs_CurveType theTarget) +{ + cf = c1; cl = c2; + Handle(Geom_Curve) c3d, newc3d[3]; + Standard_Integer i, imin = -1; + c3d = theC3d; + if (c3d.IsNull()) return newc3d[imin]; + gp_Pnt P1 = c3d->Value(c1); + gp_Pnt P2 = c3d->Value(c2); + gp_Pnt P3 = c3d->Value(c1 + (c2 - c1) / 2); + Standard_Real d[3] = { RealLast(), RealLast(), RealLast() }; + Standard_Real fp[3], lp[3]; + + if (c3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) { + Handle(Geom_TrimmedCurve) aTc = Handle(Geom_TrimmedCurve)::DownCast(c3d); + c3d = aTc->BasisCurve(); + } + + if (theConvType == GeomConvert_Target) + { + theGap = RealLast(); + if (theTarget == GeomAbs_Line) + { + newc3d[0] = ComputeLine(c3d, tolerance, c1, c2, fp[0], lp[0], theGap); + cf = fp[0]; + cl = lp[0]; + return newc3d[0]; + } + if (theTarget == GeomAbs_Circle) + { + newc3d[1] = ComputeCircle(c3d, tolerance, c1, c2, fp[1], lp[1], theGap); + cf = fp[1]; + cl = lp[1]; + return newc3d[1]; + } + if (theTarget == GeomAbs_Ellipse) + { + newc3d[2] = ComputeEllipse(c3d, tolerance, c1, c2, fp[2], lp[2], theGap); + cf = fp[2]; + cl = lp[2]; + return newc3d[2]; + } + } + // + if (theConvType == GeomConvert_Simplest) + { + theGap = RealLast(); + newc3d[0] = ComputeLine(c3d, tolerance, c1, c2, fp[0], lp[0], theGap); + if (!newc3d[0].IsNull()) + { + cf = fp[0]; + cl = lp[0]; + return newc3d[0]; + } + theGap = RealLast(); + newc3d[1] = ComputeCircle(c3d, tolerance, c1, c2, fp[1], lp[1], theGap); + if (!newc3d[1].IsNull()) + { + cf = fp[1]; + cl = lp[1]; + return newc3d[1]; + } + theGap = RealLast(); + newc3d[2] = ComputeEllipse(c3d, tolerance, c1, c2, fp[2], lp[2], theGap); + if (!newc3d[2].IsNull()) + { + cf = fp[2]; + cl = lp[2]; + return newc3d[2]; + } + // Conversion failed, returns null curve + return newc3d[0]; + } + + // theConvType == GeomConvert_MinGap + // recognition in case of small curve + imin = -1; + if((P1.Distance(P2) < 2*tolerance) && (P1.Distance(P3) < 2*tolerance)) { + newc3d[1] = ComputeCircle(c3d, tolerance, c1, c2, fp[1], lp[1], d[1]); + newc3d[0] = ComputeLine(c3d, tolerance, c1, c2, fp[0], lp[0], d[0]); + imin = 1; + if (newc3d[1].IsNull() || d[0] < d[1]) + { + imin = 0; + } + } + else { + d[0] = RealLast(); + newc3d[0] = ComputeLine (c3d,tolerance,c1,c2,fp[0],lp[0],d[0]); + Standard_Real tol = Min(tolerance, d[0]); + if (!Precision::IsInfinite(c1) && !Precision::IsInfinite(c2)) + { + d[1] = RealLast(); + newc3d[1] = ComputeCircle(c3d, tol, c1, c2, fp[1], lp[1], d[1]); + tol = Min(tol, d[1]); + d[2] = RealLast(); + newc3d[2] = ComputeEllipse(c3d, tol, c1, c2, fp[2], lp[2], d[2]); + } + Standard_Real dd = RealLast(); + for (i = 0; i < 3; ++i) + { + if (newc3d[i].IsNull()) continue; + if (d[i] < dd) + { + dd = d[i]; + imin = i; + } + } + } + + if (imin >= 0) + { + cf = fp[imin]; + cl = lp[imin]; + theGap = d[imin]; + return newc3d[imin]; + } + else + { + cf = c1; + cl = c2; + theGap = -1.; + return newc3d[0]; // must be null curve; + } +} diff --git a/src/GeomConvert/GeomConvert_CurveToAnaCurve.hxx b/src/GeomConvert/GeomConvert_CurveToAnaCurve.hxx new file mode 100644 index 0000000000..5bffa86d57 --- /dev/null +++ b/src/GeomConvert/GeomConvert_CurveToAnaCurve.hxx @@ -0,0 +1,133 @@ +// Created: 2001-05-21 +// +// Copyright (c) 2001-2013 OPEN CASCADE SAS +// +// This file is part of commercial software by OPEN CASCADE SAS, +// furnished in accordance with the terms and conditions of the contract +// and with the inclusion of this copyright notice. +// This file or any part thereof may not be provided or otherwise +// made available to any third party. +// +// No ownership title to the software is transferred hereby. +// +// OPEN CASCADE SAS makes no representation or warranties with respect to the +// performance of this software, and specifically disclaims any responsibility +// for any damages, special or consequential, connected with its use. + +#ifndef _GeomConvert_CurveToAnaCurve_HeaderFile +#define _GeomConvert_CurveToAnaCurve_HeaderFile + +#include +#include +#include +#include +#include +#include +#include +#include + +class Geom_Curve; +class Geom_Line; +class gp_Lin; +class gp_Pnt; +class gp_Circ; + + + +class GeomConvert_CurveToAnaCurve +{ +public: + + DEFINE_STANDARD_ALLOC + + Standard_EXPORT GeomConvert_CurveToAnaCurve(); + + Standard_EXPORT GeomConvert_CurveToAnaCurve(const Handle(Geom_Curve)& C); + + Standard_EXPORT void Init (const Handle(Geom_Curve)& C); + + //! Converts me to analytical if possible with given + //! tolerance. The new first and last parameters are + //! returned to newF, newL + Standard_EXPORT Standard_Boolean ConvertToAnalytical (const Standard_Real theTol, Handle(Geom_Curve)& theResultCurve, const Standard_Real F, const Standard_Real L, Standard_Real& newF, Standard_Real& newL); + + Standard_EXPORT static Handle(Geom_Curve) ComputeCurve (const Handle(Geom_Curve)& curve, const Standard_Real tolerance, + const Standard_Real c1, const Standard_Real c2, Standard_Real& cf, Standard_Real& cl, + Standard_Real& theGap, const GeomConvert_ConvType theCurvType = GeomConvert_MinGap, const GeomAbs_CurveType theTarget = GeomAbs_Line); + + //! Tries to convert the given curve to circle with given + //! tolerance. Returns NULL curve if conversion is + //! not possible. + Standard_EXPORT static Handle(Geom_Curve) ComputeCircle (const Handle(Geom_Curve)& curve, const Standard_Real tolerance, const Standard_Real c1, const Standard_Real c2, Standard_Real& cf, Standard_Real& cl, Standard_Real& Deviation); + + //! Tries to convert the given curve to ellipse with given + //! tolerance. Returns NULL curve if conversion is + //! not possible. + Standard_EXPORT static Handle(Geom_Curve) ComputeEllipse (const Handle(Geom_Curve)& curve, const Standard_Real tolerance, const Standard_Real c1, const Standard_Real c2, Standard_Real& cf, Standard_Real& cl, Standard_Real& Deviation); + + //! Tries to convert the given curve to line with given + //! tolerance. Returns NULL curve if conversion is + //! not possible. + Standard_EXPORT static Handle(Geom_Line) ComputeLine (const Handle(Geom_Curve)& curve, const Standard_Real tolerance, const Standard_Real c1, const Standard_Real c2, Standard_Real& cf, Standard_Real& cl, Standard_Real& Deviation); + + //! Returns true if the set of points is linear with given + //! tolerance + Standard_EXPORT static Standard_Boolean IsLinear (const TColgp_Array1OfPnt& aPoints, const Standard_Real tolerance, Standard_Real& Deviation); + + //! Creates line on two points. + //! Resulting parameters returned + Standard_EXPORT static gp_Lin GetLine(const gp_Pnt& P1, const gp_Pnt& P2, Standard_Real& cf, Standard_Real& cl); + + //! Creates circle on points. Returns true if OK. + Standard_EXPORT static Standard_Boolean GetCircle(gp_Circ& Circ, const gp_Pnt& P0, const gp_Pnt& P1, const gp_Pnt& P2); + + //! Returns maximal deviation of converted surface from the original + //! one computed by last call to ConvertToAnalytical + Standard_Real Gap() const + { + return myGap; + } + + //! Returns conversion type + GeomConvert_ConvType GetConvType() const + { + return myConvType; + } + + //! Sets type of convertion + void SetConvType(const GeomConvert_ConvType theConvType) + { + myConvType = theConvType; + } + + //! Returns target curve type + GeomAbs_CurveType GetTarget() const + { + return myTarget; + } + + //! Sets target curve type + void SetTarget(const GeomAbs_CurveType theTarget) + { + myTarget = theTarget; + } + + +protected: + + + + + +private: + + Handle(Geom_Curve) myCurve; + Standard_Real myGap; + GeomConvert_ConvType myConvType; + GeomAbs_CurveType myTarget; + +}; + + + +#endif // _GeomConvert_CurveToAnaCurve_HeaderFile diff --git a/src/GeomConvert/GeomConvert_SurfToAnaSurf.cxx b/src/GeomConvert/GeomConvert_SurfToAnaSurf.cxx new file mode 100644 index 0000000000..37ed349ebb --- /dev/null +++ b/src/GeomConvert/GeomConvert_SurfToAnaSurf.cxx @@ -0,0 +1,764 @@ +// Created: 1998-06-03 +// +// Copyright (c) 1999-2013 OPEN CASCADE SAS +// +// This file is part of commercial software by OPEN CASCADE SAS, +// furnished in accordance with the terms and conditions of the contract +// and with the inclusion of this copyright notice. +// This file or any part thereof may not be provided or otherwise +// made available to any third party. +// +// No ownership title to the software is transferred hereby. +// +// OPEN CASCADE SAS makes no representation or warranties with respect to the +// performance of this software, and specifically disclaims any responsibility +// for any damages, special or consequential, connected with its use. + +//abv 06.01.99 fix of misprint +//:p6 abv 26.02.99: make ConvertToPeriodic() return Null if nothing done + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +//static method for checking surface of revolution +//To avoid two-parts cone-like surface +static void CheckVTrimForRevSurf(const Handle(Geom_SurfaceOfRevolution)& aRevSurf, + Standard_Real& V1, Standard_Real& V2) +{ + const Handle(Geom_Curve)& aBC = aRevSurf->BasisCurve(); + Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(aBC); + if (aLine.IsNull()) + return; + const gp_Ax1& anAxis = aRevSurf->Axis(); + + gp_Lin anALin(anAxis); + Extrema_ExtElC anExtLL(aLine->Lin(), anALin, Precision::Angular()); + if (!anExtLL.IsDone() || anExtLL.IsParallel()) + return; + Standard_Integer aNbExt = anExtLL.NbExt(); + if (aNbExt == 0) + return; + + Standard_Integer i; + Standard_Integer imin = 0; + for (i = 1; i <= aNbExt; ++i) + { + if (anExtLL.SquareDistance(i) < Precision::SquareConfusion()) + { + imin = i; + break; + } + } + if (imin == 0) + return; + + Extrema_POnCurv aP1, aP2; + anExtLL.Points(imin, aP1, aP2); + Standard_Real aVExt = aP1.Parameter(); + if (aVExt <= V1 || aVExt >= V2) + return; + + if (aVExt - V1 > V2 - aVExt) + { + V2 = aVExt; + } + else + { + V1 = aVExt; + } + +} + +// static method to try create toroidal surface. +// In case = Standard_True try to use V isoline radius as minor radaius. +static Handle(Geom_Surface) TryTorusSphere(const Handle(Geom_Surface)& theSurf, + const Handle(Geom_Circle)& circle, + const Handle(Geom_Circle)& otherCircle, + const Standard_Real Param1, + const Standard_Real Param2, + const Standard_Real aParam1ToCrv, + const Standard_Real aParam2ToCrv, + const Standard_Real toler, + const Standard_Boolean isTryUMajor) +{ + Handle(Geom_Surface) newSurface; + Standard_Real cf, cl; + Handle(Geom_Curve) IsoCrv1; + Handle(Geom_Curve) IsoCrv2; + Standard_Real aGap1, aGap2; + // initial radius + Standard_Real R = circle->Circ().Radius(); + // iso lines + + if (isTryUMajor) + { + IsoCrv1 = theSurf->VIso(Param1 + ((Param2 - Param1) / 3.)); + IsoCrv2 = theSurf->VIso(Param1 + ((Param2 - Param1)*2. / 3)); + } + else + { + IsoCrv1 = theSurf->UIso(Param1 + ((Param2 - Param1) / 3.)); + IsoCrv2 = theSurf->UIso(Param1 + ((Param2 - Param1)*2. / 3)); + } + + Handle(Geom_Curve) Crv1 = GeomConvert_CurveToAnaCurve::ComputeCurve(IsoCrv1, toler, aParam1ToCrv, aParam2ToCrv, cf, cl, + aGap1); + Handle(Geom_Curve) Crv2 = GeomConvert_CurveToAnaCurve::ComputeCurve(IsoCrv2, toler, aParam1ToCrv, aParam2ToCrv, cf, cl, + aGap2); + if (Crv1.IsNull() || Crv2.IsNull() || + !Crv1->IsKind(STANDARD_TYPE(Geom_Circle)) || + !Crv2->IsKind(STANDARD_TYPE(Geom_Circle))) + return newSurface; + + Handle(Geom_Circle) aCircle1 = Handle(Geom_Circle)::DownCast(Crv1); + Handle(Geom_Circle) aCircle2 = Handle(Geom_Circle)::DownCast(Crv2); + Standard_Real R1 = aCircle1->Circ().Radius(); + Standard_Real R2 = aCircle2->Circ().Radius(); + + // check radiuses + if ((Abs(R - R1) > toler) || (Abs(R - R2) > toler)) + return newSurface; + + // get centers of the major radius + gp_Pnt aPnt1, aPnt2, aPnt3; + aPnt1 = circle->Circ().Location(); + aPnt2 = aCircle1->Circ().Location(); + aPnt3 = aCircle2->Circ().Location(); + + //Standard_Real eps = 1.e-09; // angular resolution + Standard_Real d0 = aPnt1.Distance(aPnt2); Standard_Real d1 = aPnt1.Distance(aPnt3); + gp_Circ circ; + + if (d0 < toler || d1 < toler) + { + // compute sphere + gp_Dir MainDir = otherCircle->Circ().Axis().Direction(); + gp_Ax3 Axes(circle->Circ().Location(), MainDir); + Handle(Geom_SphericalSurface) anObject = + new Geom_SphericalSurface(Axes, R); + if (!anObject.IsNull()) + newSurface = anObject; + + return newSurface; + } + + if (!GeomConvert_CurveToAnaCurve::GetCircle(circ, aPnt1, aPnt2, aPnt3) /*, d0, d1, eps)*/) + return newSurface; + + Standard_Real aMajorR = circ.Radius(); + gp_Pnt aCenter = circ.Location(); + gp_Dir aDir((aPnt1.XYZ() - aCenter.XYZ()) ^ (aPnt3.XYZ() - aCenter.XYZ())); + gp_Ax3 anAx3(aCenter, aDir); + newSurface = new Geom_ToroidalSurface(anAx3, aMajorR, R); + return newSurface; +} + +static Standard_Real ComputeGap(const Handle(Geom_Surface)& theSurf, + const Standard_Real theU1, const Standard_Real theU2, const Standard_Real theV1, const Standard_Real theV2, + const Handle(Geom_Surface) theNewSurf, const Standard_Real theTol = RealLast()) +{ + GeomAdaptor_Surface aGAS(theNewSurf); + GeomAbs_SurfaceType aSType = aGAS.GetType(); + gp_Pln aPln; + gp_Cylinder aCyl; + gp_Cone aCon; + gp_Sphere aSphere; + gp_Torus aTor; + switch (aSType) + { + case GeomAbs_Plane: + aPln = aGAS.Plane(); + break; + case GeomAbs_Cylinder: + aCyl = aGAS.Cylinder(); + break; + case GeomAbs_Cone: + aCon = aGAS.Cone(); + break; + case GeomAbs_Sphere: + aSphere = aGAS.Sphere(); + break; + case GeomAbs_Torus: + aTor = aGAS.Torus(); + break; + default: + return Precision::Infinite(); + break; + } + + Standard_Real aGap = 0.; + Standard_Boolean onSurface = Standard_True; + + Standard_Real S, T; + gp_Pnt P3d, P3d2; + + const Standard_Integer NP = 21; + Standard_Real DU, DV; + Standard_Integer j, i; + DU = (theU2 - theU1) / (NP - 1); + DV = (theV2 - theV1) / (NP - 1); + Standard_Real DU2 = DU / 2., DV2 = DV / 2.; + for (j = 1; (j < NP) && onSurface; j++) { + Standard_Real V = theV1 + DV*(j - 1) + DV2; + for (i = 1; i <= NP; i++) { + Standard_Real U = theU1 + DU*(i - 1) + DU2; + theSurf->D0(U, V, P3d); + + switch (aSType) { + + case GeomAbs_Plane: + { + ElSLib::Parameters(aPln, P3d, S, T); + P3d2 = ElSLib::Value(S, T, aPln); + break; + } + case GeomAbs_Cylinder: + { + ElSLib::Parameters(aCyl, P3d, S, T); + P3d2 = ElSLib::Value(S, T, aCyl); + break; + } + case GeomAbs_Cone: + { + ElSLib::Parameters(aCon, P3d, S, T); + P3d2 = ElSLib::Value(S, T, aCon); + break; + } + case GeomAbs_Sphere: + { + ElSLib::Parameters(aSphere, P3d, S, T); + P3d2 = ElSLib::Value(S, T, aSphere); + break; + } + case GeomAbs_Torus: + { + ElSLib::Parameters(aTor, P3d, S, T); + P3d2 = ElSLib::Value(S, T, aTor); + break; + } + default: + S = 0.; T = 0.; + theNewSurf->D0(S, T, P3d2); + break; + } + + Standard_Real dis = P3d.Distance(P3d2); + if (dis > aGap) aGap = dis; + + if (aGap > theTol) { + onSurface = Standard_False; + break; + } + } + } + return aGap; +} +//======================================================================= +//function : GeomConvert_SurfToAnaSurf +//purpose : +//======================================================================= +GeomConvert_SurfToAnaSurf::GeomConvert_SurfToAnaSurf() : + myGap (-1.), + myConvType (GeomConvert_Simplest), + myTarget (GeomAbs_Plane) +{ +} + +//======================================================================= +//function : GeomConvert_SurfToAnaSurf +//purpose : +//======================================================================= + +GeomConvert_SurfToAnaSurf::GeomConvert_SurfToAnaSurf (const Handle(Geom_Surface)& S) : + myGap(-1.), + myConvType(GeomConvert_Simplest), + myTarget(GeomAbs_Plane) +{ + Init ( S ); +} + +//======================================================================= +//function : Init +//purpose : +//======================================================================= + +void GeomConvert_SurfToAnaSurf::Init (const Handle(Geom_Surface)& S) +{ + mySurf = S; +} + +//======================================================================= +//function : ConvertToAnalytical +//purpose : +//======================================================================= + +Handle(Geom_Surface) GeomConvert_SurfToAnaSurf::ConvertToAnalytical (const Standard_Real InitialToler) +{ + Standard_Real U1, U2, V1, V2; + mySurf->Bounds(U1, U2, V1, V2); + if(Precision::IsInfinite(U1) && Precision::IsInfinite(U2) ) { + U1 = -1.; + U2 = 1.; + } + if(Precision::IsInfinite(V1) && Precision::IsInfinite(V2) ) { + V1 = -1.; + V2 = 1.; + Handle(Geom_SurfaceOfRevolution) aRevSurf = Handle(Geom_SurfaceOfRevolution)::DownCast(mySurf); + if (!aRevSurf.IsNull()) + { + CheckVTrimForRevSurf(aRevSurf, V1, V2); + } + } + return ConvertToAnalytical(InitialToler, U1, U2, V1, V2); +} + +//======================================================================= +//function : ConvertToAnalytical +//purpose : +//======================================================================= +Handle(Geom_Surface) GeomConvert_SurfToAnaSurf::ConvertToAnalytical(const Standard_Real InitialToler, + const Standard_Real Umin, const Standard_Real Umax, + const Standard_Real Vmin, const Standard_Real Vmax) +{ + // + GeomAdaptor_Surface aGAS(mySurf); + GeomAbs_SurfaceType aSType = aGAS.GetType(); + switch (aSType) + { + case GeomAbs_Plane: + { + myGap = 0.; + return new Geom_Plane(aGAS.Plane()); + } + case GeomAbs_Cylinder: + { + myGap = 0.; + return new Geom_CylindricalSurface(aGAS.Cylinder()); + } + case GeomAbs_Cone: + { + myGap = 0.; + return new Geom_ConicalSurface(aGAS.Cone()); + } + case GeomAbs_Sphere: + { + myGap = 0.; + return new Geom_SphericalSurface(aGAS.Sphere()); + } + case GeomAbs_Torus: + { + myGap = 0.; + return new Geom_ToroidalSurface(aGAS.Torus()); + } + default: + break; + } + // + Standard_Real toler = InitialToler; + Handle(Geom_Surface) newSurf[5]; + Standard_Real dd[5] = { -1., -1., -1., -1., -1 }; + + //Check boundaries + Standard_Real U1, U2, V1, V2; + mySurf->Bounds(U1, U2, V1, V2); + Standard_Boolean aDoSegment = Standard_False; + Standard_Real aTolBnd = Precision::PConfusion(); + Standard_Integer isurf = 0; + if (Umin < U1 || Umax > U2 || Vmin < V1 || Vmax > V2) + { + return newSurf[isurf]; + } + else + { + if (Umin - U1 > aTolBnd) + { + U1 = Umin; + aDoSegment = Standard_True; + } + if (U2 - Umax > aTolBnd) + { + U2 = Umax; + aDoSegment = Standard_True; + } + if (Vmin - V1 > aTolBnd) + { + V1 = Vmin; + aDoSegment = Standard_True; + } + if (V2 - Vmax > aTolBnd) + { + V2 = Vmax; + aDoSegment = Standard_True; + } + + } + + Standard_Boolean IsBz = aSType == GeomAbs_BezierSurface; + Standard_Boolean IsBs = aSType == GeomAbs_BSplineSurface; + + Handle(Geom_Surface) aTempS = mySurf; + if (IsBs) + { + Handle(Geom_BSplineSurface) aBs = Handle(Geom_BSplineSurface)::DownCast(mySurf->Copy()); + if (aDoSegment) + { + aBs->Segment(U1, U2, V1, V2); + } + aTempS = aBs; + } + else if (IsBz) + { + Handle(Geom_BezierSurface) aBz = Handle(Geom_BezierSurface)::DownCast(mySurf->Copy()); + if (aDoSegment) + { + aBz->Segment(U1, U2, V1, V2); + } + aTempS = aBz; + } + // check the planarity first + if (!IsBs && !IsBz) + { + aTempS = new Geom_RectangularTrimmedSurface(aTempS, U1, U2, V1, V2); + } + isurf = 0; // set plane + GeomLib_IsPlanarSurface GeomIsPlanar(aTempS, toler); + if (GeomIsPlanar.IsPlanar()) { + gp_Pln newPln = GeomIsPlanar.Plan(); + newSurf[isurf] = new Geom_Plane(newPln); + dd[isurf] = ComputeGap(aTempS, U1, U2, V1, V2, newSurf[isurf]); + if ( myConvType == GeomConvert_Simplest || + (myConvType == GeomConvert_Target && myTarget == GeomAbs_Plane)) + { + myGap = dd[isurf]; + return newSurf[isurf]; + } + } + else + { + if (myConvType == GeomConvert_Target && myTarget == GeomAbs_Plane) + { + myGap = dd[isurf]; + return newSurf[isurf]; + } + } + + Standard_Real diagonal = mySurf->Value(U1, V1).Distance(mySurf->Value((U1 + U2), (V1 + V2) / 2)); + Standard_Real twist = 1000; + if (toler > diagonal / twist) toler = diagonal / twist; + + isurf = 1; // set cylinder + Standard_Boolean aCylinderConus = Standard_False; + Standard_Boolean aToroidSphere = Standard_False; + + //convert middle uiso and viso to canonical representation + Standard_Real VMid = 0.5*(V1 + V2); + Standard_Real UMid = 0.5*(U1 + U2); + Handle(Geom_Surface) TrSurf = aTempS; + if(!aDoSegment) + TrSurf = new Geom_RectangularTrimmedSurface(aTempS, U1, U2, V1, V2); + + Handle(Geom_Curve) UIso = TrSurf->UIso(UMid); + Handle(Geom_Curve) VIso = TrSurf->VIso(VMid); + + Standard_Real cuf, cul, cvf, cvl, aGap1, aGap2; + Handle(Geom_Curve) umidiso = GeomConvert_CurveToAnaCurve::ComputeCurve(UIso, toler, V1, V2, cuf, cul, aGap1); + Handle(Geom_Curve) vmidiso = GeomConvert_CurveToAnaCurve::ComputeCurve(VIso, toler, U1, U2, cvf, cvl, aGap2); + if (umidiso.IsNull() || vmidiso.IsNull()) { + return newSurf[isurf]; + } + + // + Standard_Boolean VCase = Standard_False; + + if (umidiso->IsKind(STANDARD_TYPE(Geom_Circle)) && vmidiso->IsKind(STANDARD_TYPE(Geom_Circle))) + { + aToroidSphere = Standard_True; + if (myConvType == GeomConvert_Target && (myTarget == GeomAbs_Cylinder || myTarget == GeomAbs_Cone)) + { + isurf = 1; + myGap = dd[isurf]; + return newSurf[isurf]; + } + isurf = 3; // set sphere + } + else if (umidiso->IsKind(STANDARD_TYPE(Geom_Line)) && vmidiso->IsKind(STANDARD_TYPE(Geom_Circle))) { + aCylinderConus = Standard_True; VCase = Standard_True; + if (myConvType == GeomConvert_Target && (myTarget == GeomAbs_Sphere || myTarget == GeomAbs_Torus)) + { + isurf = 3; + myGap = dd[isurf]; + return newSurf[isurf]; + } + isurf = 1;// set cylinder + } + else if (umidiso->IsKind(STANDARD_TYPE(Geom_Circle)) && vmidiso->IsKind(STANDARD_TYPE(Geom_Line))) + { + aCylinderConus = Standard_True; + if (myConvType == GeomConvert_Target && (myTarget == GeomAbs_Sphere || myTarget == GeomAbs_Torus)) + { + isurf = 3; + myGap = dd[isurf]; + return newSurf[isurf]; + } + isurf = 1;// set cylinder + } + Standard_Real cl = 0.0; + //case of torus-sphere + if (aToroidSphere) { + + isurf = 3; // Set spherical surface + + Handle(Geom_Circle) Ucircle = Handle(Geom_Circle)::DownCast(umidiso); + Handle(Geom_Circle) Vcircle = Handle(Geom_Circle)::DownCast(vmidiso); + // torus + // try when V isolines is with same radius + Handle(Geom_Surface) anObject = + TryTorusSphere(mySurf, Vcircle, Ucircle, V1, V2, U1, U2, toler, Standard_True); + if (anObject.IsNull()) // try when U isolines is with same radius + anObject = TryTorusSphere(mySurf, Ucircle, Vcircle, U1, U2, V1, V2, toler, Standard_False); + + if (!anObject.IsNull()) + { + if (anObject->IsKind(STANDARD_TYPE(Geom_ToroidalSurface))) + { + isurf = 4; // set torus + } + newSurf[isurf] = anObject; + } + else + { + myGap = dd[isurf]; + } + } + //case of cone - cylinder + else if (aCylinderConus) { + Standard_Real param1, param2, cf1, cf2; + Handle(Geom_Curve) firstiso, lastiso; + Handle(Geom_Circle) firstisocirc, lastisocirc, midisocirc; + gp_Dir isoline; + if (VCase) { + param1 = U1; param2 = U2; + firstiso = TrSurf->VIso(V1); + lastiso = TrSurf->VIso(V2); + midisocirc = Handle(Geom_Circle)::DownCast(vmidiso); + isoline = Handle(Geom_Line)::DownCast(umidiso)->Lin().Direction(); + } + else { + param1 = V1; param2 = V2; + firstiso = TrSurf->UIso(U1); + lastiso = TrSurf->UIso(U2); + midisocirc = Handle(Geom_Circle)::DownCast(umidiso); + isoline = Handle(Geom_Line)::DownCast(vmidiso)->Lin().Direction(); + } + firstisocirc = Handle(Geom_Circle)::DownCast(GeomConvert_CurveToAnaCurve::ComputeCurve(firstiso, toler, param1, param2, cf1, cl, aGap1)); + lastisocirc = Handle(Geom_Circle)::DownCast(GeomConvert_CurveToAnaCurve::ComputeCurve(lastiso, toler, param1, param2, cf2, cl, aGap2)); + if (!firstisocirc.IsNull() || !lastisocirc.IsNull()) { + Standard_Real R1, R2, R3; + gp_Pnt P1, P2, P3; + if (!firstisocirc.IsNull()) { + R1 = firstisocirc->Circ().Radius(); + P1 = firstisocirc->Circ().Location(); + } + else { + R1 = 0; + P1 = firstiso->Value((firstiso->LastParameter() - firstiso->FirstParameter()) / 2); + } + R2 = midisocirc->Circ().Radius(); + P2 = midisocirc->Circ().Location(); + if (!lastisocirc.IsNull()) { + R3 = lastisocirc->Circ().Radius(); + P3 = lastisocirc->Circ().Location(); + } + else { + R3 = 0; + P3 = lastiso->Value((lastiso->LastParameter() - lastiso->FirstParameter()) / 2); + } + //cylinder + if (((Abs(R2 - R1)) < toler) && ((Abs(R3 - R1)) < toler) && + ((Abs(R3 - R2)) < toler)) { + isurf = 1; + gp_Ax3 Axes(P1, gp_Dir(gp_Vec(P1, P3))); + Handle(Geom_CylindricalSurface) anObject = + new Geom_CylindricalSurface(Axes, R1); + if (myConvType == GeomConvert_Target && myTarget != GeomAbs_Cylinder) + { + return newSurf[isurf]; + } + newSurf[isurf] = anObject; + } + //cone + else if ((((Abs(R1)) > (Abs(R2))) && ((Abs(R2)) > (Abs(R3)))) || + (((Abs(R3)) > (Abs(R2))) && ((Abs(R2)) > (Abs(R1))))) { + Standard_Real radius; + gp_Ax3 Axes; + Standard_Real semiangle = + gp_Vec(isoline).Angle(gp_Vec(P3, P1)); + if (semiangle>M_PI / 2) semiangle = M_PI - semiangle; + if (R1 > R3) { + radius = R3; + Axes = gp_Ax3(P3, gp_Dir(gp_Vec(P3, P1))); + } + else { + radius = R1; + Axes = gp_Ax3(P1, gp_Dir(gp_Vec(P1, P3))); + } + Handle(Geom_ConicalSurface) anObject = + new Geom_ConicalSurface(Axes, semiangle, radius); + isurf = 2; + if (myConvType == GeomConvert_Target && myTarget != GeomAbs_Cone) + { + return newSurf[isurf]; + } + newSurf[isurf] = anObject; + } + } + else + { + myGap = dd[isurf]; + } + } + // + //--------------------------------------------------------------------- + // verification + //--------------------------------------------------------------------- + GeomAbs_SurfaceType aSTypes[5] = { GeomAbs_Plane, GeomAbs_Cylinder, + GeomAbs_Cone, GeomAbs_Sphere, GeomAbs_Torus }; + Standard_Integer imin = -1; + Standard_Real aDmin = RealLast(); + for (isurf = 0; isurf < 5; ++isurf) + { + if (newSurf[isurf].IsNull()) + continue; + dd[isurf] = ComputeGap(TrSurf, U1, U2, V1, V2, newSurf[isurf], toler); + if (dd[isurf] <= toler) + { + if (myConvType == GeomConvert_Simplest || (myConvType == GeomConvert_Target && myTarget == aSTypes[isurf])) + { + myGap = dd[isurf]; + return newSurf[isurf]; + } + else if (myConvType == GeomConvert_MinGap) + { + if (dd[isurf] < aDmin) + { + aDmin = dd[isurf]; + imin = isurf; + } + } + } + } + // + if (imin >= 0) + { + myGap = dd[imin]; + return newSurf[imin]; + } + + return newSurf[0]; +} + +//======================================================================= +//function : IsSame +//purpose : +//======================================================================= + +Standard_Boolean GeomConvert_SurfToAnaSurf::IsSame(const Handle(Geom_Surface)& S1, + const Handle(Geom_Surface)& S2, + const Standard_Real tol) +{ + // only elementary surfaces are processed + if(!S1->IsKind(STANDARD_TYPE(Geom_ElementarySurface)) || + !S2->IsKind(STANDARD_TYPE(Geom_ElementarySurface))) + return Standard_False; + + Handle(GeomAdaptor_Surface) anAdaptor1 = new GeomAdaptor_Surface(S1); + Handle(GeomAdaptor_Surface) anAdaptor2 = new GeomAdaptor_Surface(S2); + + GeomAbs_SurfaceType aST1 = anAdaptor1->GetType(); + GeomAbs_SurfaceType aST2 = anAdaptor2->GetType(); + + if (aST1 != aST2) + { + return Standard_False; + } + + IntAna_QuadQuadGeo interii; + if (aST1 == GeomAbs_Plane) + { + interii.Perform(anAdaptor1->Plane(), anAdaptor2->Plane(), tol, tol); + } + else if (aST1 == GeomAbs_Cylinder) + { + interii.Perform(anAdaptor1->Cylinder(), anAdaptor2->Cylinder(), tol); + } + else if (aST1 == GeomAbs_Cone) + { + interii.Perform(anAdaptor1->Cone(), anAdaptor2->Cone(), tol); + } + else if (aST1 == GeomAbs_Sphere) + { + interii.Perform(anAdaptor1->Sphere(), anAdaptor2->Sphere(), tol); + } + else if (aST1 == GeomAbs_Torus) + { + interii.Perform(anAdaptor1->Torus(), anAdaptor2->Torus(), tol); + } + + if (!interii.IsDone()) + return Standard_False; + + IntAna_ResultType aTypeRes = interii.TypeInter(); + + return aTypeRes == IntAna_Same; + +} + + +//======================================================================= +//function : IsCanonical +//purpose : +//======================================================================= + +Standard_Boolean GeomConvert_SurfToAnaSurf::IsCanonical(const Handle(Geom_Surface)& S) +{ + if(S.IsNull()) + return Standard_False; + + if ( S->IsKind(STANDARD_TYPE(Geom_Plane)) + || S->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)) + || S->IsKind(STANDARD_TYPE(Geom_ConicalSurface)) + || S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) + || S->IsKind(STANDARD_TYPE(Geom_ToroidalSurface))) + return Standard_True; + + return Standard_False; +} + diff --git a/src/GeomConvert/GeomConvert_SurfToAnaSurf.hxx b/src/GeomConvert/GeomConvert_SurfToAnaSurf.hxx new file mode 100644 index 0000000000..016956590e --- /dev/null +++ b/src/GeomConvert/GeomConvert_SurfToAnaSurf.hxx @@ -0,0 +1,91 @@ +// Created: 1998-06-03 +// +// Copyright (c) 1999-2013 OPEN CASCADE SAS +// +// This file is part of commercial software by OPEN CASCADE SAS, +// furnished in accordance with the terms and conditions of the contract +// and with the inclusion of this copyright notice. +// This file or any part thereof may not be provided or otherwise +// made available to any third party. +// +// No ownership title to the software is transferred hereby. +// +// OPEN CASCADE SAS makes no representation or warranties with respect to the +// performance of this software, and specifically disclaims any responsibility +// for any damages, special or consequential, connected with its use. + +#ifndef _GeomConvert_SurfToAnaSurf_HeaderFile +#define _GeomConvert_SurfToAnaSurf_HeaderFile + +#include +#include +#include + +#include +#include +#include +#include +class Geom_Surface; + + +//! Converts a surface to the analitical form with given +//! precision. Conversion is done only the surface is bspline +//! of bezier and this can be approximed by some analytical +//! surface with that precision. +class GeomConvert_SurfToAnaSurf +{ +public: + + DEFINE_STANDARD_ALLOC + + Standard_EXPORT GeomConvert_SurfToAnaSurf(); + + Standard_EXPORT GeomConvert_SurfToAnaSurf(const Handle(Geom_Surface)& S); + + Standard_EXPORT void Init (const Handle(Geom_Surface)& S); + + void SetConvType(const GeomConvert_ConvType theConvType = GeomConvert_Simplest) + { + myConvType = theConvType; + } + void SetTarget(const GeomAbs_SurfaceType theSurfType = GeomAbs_Plane) + { + myTarget = theSurfType; + } + + //! Returns maximal deviation of converted surface from the original + //! one computed by last call to ConvertToAnalytical + Standard_Real Gap() const + { + return myGap; + } + + //! Tries to convert the Surface to an Analytic form + //! Returns the result + //! In case of failure, returns a Null Handle + //! + Standard_EXPORT Handle(Geom_Surface) ConvertToAnalytical (const Standard_Real InitialToler); + Standard_EXPORT Handle(Geom_Surface) ConvertToAnalytical (const Standard_Real InitialToler, + const Standard_Real Umin, const Standard_Real Umax, + const Standard_Real Vmin, const Standard_Real Vmax); + + //! Returns true if surfaces is same with the given tolerance + Standard_EXPORT static Standard_Boolean IsSame (const Handle(Geom_Surface)& S1, const Handle(Geom_Surface)& S2, const Standard_Real tol); + + //! Returns true, if surface is canonical + Standard_EXPORT static Standard_Boolean IsCanonical (const Handle(Geom_Surface)& S); + + +protected: + +private: + + Handle(Geom_Surface) mySurf; + Standard_Real myGap; + GeomConvert_ConvType myConvType; + GeomAbs_SurfaceType myTarget; + +}; + + +#endif // _GeomConvert_SurfToAnaSurf_HeaderFile diff --git a/src/GeomliteTest/GeomliteTest_SurfaceCommands.cxx b/src/GeomliteTest/GeomliteTest_SurfaceCommands.cxx index 4becc6f98a..37bf62d118 100644 --- a/src/GeomliteTest/GeomliteTest_SurfaceCommands.cxx +++ b/src/GeomliteTest/GeomliteTest_SurfaceCommands.cxx @@ -59,7 +59,9 @@ #include #include #include - +#include +#include +#include #include #include @@ -522,6 +524,111 @@ static Standard_Integer converting(Draw_Interpretor& , Standard_Integer n, const return 0; } +//======================================================================= +//function : converting to canonical +//purpose : +//======================================================================= + +static Standard_Integer tocanon(Draw_Interpretor& di, Standard_Integer n, const char ** a) +{ + if (n < 3) return 1; + + GeomConvert_ConvType aConvType = GeomConvert_Simplest; + GeomAbs_CurveType aCurv = GeomAbs_Line; + GeomAbs_SurfaceType aSurf = GeomAbs_Plane; + if (n > 4) + { + if (strcmp(a[4], "sim") == 0) { + aConvType = GeomConvert_Simplest; + } + else if (strcmp(a[4], "gap") == 0) { + aConvType = GeomConvert_MinGap; + } + else if (strcmp(a[4], "lin") == 0) { + aConvType = GeomConvert_Target; + aCurv = GeomAbs_Line; + } + else if (strcmp(a[4], "cir") == 0) { + aConvType = GeomConvert_Target; + aCurv = GeomAbs_Circle; + } + else if (strcmp(a[4], "ell") == 0) { + aConvType = GeomConvert_Target; + aCurv = GeomAbs_Ellipse; + } + else if (strcmp(a[4], "pln") == 0) { + aConvType = GeomConvert_Target; + aSurf = GeomAbs_Plane; + } + else if (strcmp(a[4], "cyl") == 0) { + aConvType = GeomConvert_Target; + aSurf = GeomAbs_Cylinder; + } + else if (strcmp(a[4], "con") == 0) { + aConvType = GeomConvert_Target; + aSurf = GeomAbs_Cone; + } + else if (strcmp(a[4], "sph") == 0) { + aConvType = GeomConvert_Target; + aSurf = GeomAbs_Sphere; + } + else if (strcmp(a[4], "tor") == 0) { + aConvType = GeomConvert_Target; + aSurf = GeomAbs_Torus; + } + } + + Standard_Real tol = Precision::Confusion(); + if (n > 3) + { + tol = Draw::Atof(a[3]); + } + + Handle(Geom_Curve) GC = DrawTrSurf::GetCurve(a[2]); + if (GC.IsNull()) { + Handle(Geom_Surface) GS = DrawTrSurf::GetSurface(a[2]); + if (GS.IsNull()) { + return 1; + } + else { + GeomConvert_SurfToAnaSurf aSurfToAna(GS); + aSurfToAna.SetConvType(aConvType); + if(aConvType == GeomConvert_Target) + aSurfToAna.SetTarget(aSurf); + Handle(Geom_Surface) anAnaSurf = aSurfToAna.ConvertToAnalytical(tol); + if (!anAnaSurf.IsNull()) + { + DrawTrSurf::Set(a[1], anAnaSurf); + Standard_Real aGap = aSurfToAna.Gap(); + di << "Gap = " << aGap << "\n"; + } + else + di << "Conversion failed" << "\n"; + } + } + else { + GeomConvert_CurveToAnaCurve aCurvToAna(GC); + aCurvToAna.SetConvType(aConvType); + if (aConvType == GeomConvert_Target) + aCurvToAna.SetTarget(aCurv); + + Handle(Geom_Curve) anAnaCurv; + Standard_Real tf = GC->FirstParameter(), tl = GC->LastParameter(), ntf, ntl; + Standard_Boolean isdone = aCurvToAna.ConvertToAnalytical(tol, anAnaCurv, tf, tl, ntf, ntl); + if (isdone) + { + anAnaCurv = new Geom_TrimmedCurve(anAnaCurv, ntf, ntl); + DrawTrSurf::Set(a[1], anAnaCurv); + Standard_Real aGap = aCurvToAna.Gap(); + di << "Gap = " << aGap << "\n"; + } + else + di << "Conversion failed" << "\n"; + } + + return 0; +} + //======================================================================= //function : tobezier @@ -1714,6 +1821,11 @@ void GeomliteTest::SurfaceCommands(Draw_Interpretor& theCommands) __FILE__, converting,g); + theCommands.Add("tocanon", + "tocanon result c3d/surf [tol [sim gap lin cir ell pln cyl con sph tor]]", + __FILE__, + tocanon, g); + theCommands.Add("tobezier", "tobezier result c2d/c3d/surf [ufirst, ulast / ufirst, ulast, vfirst, vlast]", __FILE__, diff --git a/src/SWDRAW/SWDRAW_ShapeAnalysis.cxx b/src/SWDRAW/SWDRAW_ShapeAnalysis.cxx index c07e225f40..26b00cbe47 100644 --- a/src/SWDRAW/SWDRAW_ShapeAnalysis.cxx +++ b/src/SWDRAW/SWDRAW_ShapeAnalysis.cxx @@ -61,6 +61,14 @@ #include #include #include +#include +#include +#include +#include +#include +#include +#include +#include #include static Standard_Integer tolerance @@ -976,6 +984,203 @@ static Standard_Integer checkedge(Draw_Interpretor& di, Standard_Integer argc, c return 0; } +//======================================================================= +// getanasurf +//======================================================================= +static Standard_Integer getanasurf(Draw_Interpretor& di, + Standard_Integer n, const char** a) + +{ + if (n < 3) { + di << "Usage: \n"; + di << "getanasurf res shape [target [tol [sample]]] \n"; + di << "target is reqired type of surface and can be: pln, cyl, con sph \n"; + di << "sample is surface of required type, which parameters are used as starting \n"; + di << "point for seaching parametrs of surface by Least Square method when input shape \n"; + di << "is edge or wire \n"; + return 1; + } + TopoDS_Shape sh = DBRep::Get(a[2]); + if (sh.IsNull()) return 1; + TopAbs_ShapeEnum aShType = sh.ShapeType(); + if (!(aShType == TopAbs_SHELL || aShType == TopAbs_FACE || aShType == TopAbs_EDGE || aShType == TopAbs_WIRE)) + { + di << "Wrong shape type, shape can be shell or face or edge or wire\n"; + return 1; + } + + GeomAbs_SurfaceType aTargets[] = { GeomAbs_Plane, GeomAbs_Cylinder, GeomAbs_Cone, GeomAbs_Sphere }; + Standard_Integer isurf = 0; + if (n > 3) + { + if (strcmp(a[3], "pln") == 0) + isurf = 0; + else if (strcmp(a[3], "cyl") == 0) + isurf = 1; + else if (strcmp(a[3], "con") == 0) + isurf = 2; + else if (strcmp(a[3], "sph") == 0) + isurf = 3; + } + + Standard_Real tol = 1.e-7; + if (n > 4) + tol = Draw::Atof(a[4]); + + // get sample target for edge and wire + GeomAdaptor_Surface aSampleSurf; + if (n > 5 && (sh.ShapeType() == TopAbs_EDGE || sh.ShapeType() == TopAbs_WIRE )) + { + Handle(Geom_Surface) aGSurf = DrawTrSurf::GetSurface(a[5]); + if (aGSurf.IsNull()) + { + di << "Sample surface is null" << "\n"; + return 1; + } + aSampleSurf.Load(aGSurf); + GeomAbs_SurfaceType aSType = aSampleSurf.GetType(); + if (aSType != aTargets[isurf]) + { + di << "Sample surface has wrong type" << "\n"; + return 1; + } + } + + ShapeAnalysis_CanonicalRecognition aCanonRec(sh); + Handle(Geom_Surface) aRes; + switch (aTargets[isurf]) + { + case GeomAbs_Plane: + { + gp_Pln aPln; + if (aSampleSurf.GetType() == GeomAbs_Plane) + aPln = aSampleSurf.Plane(); + if (aCanonRec.IsPlane(tol, aPln)) + aRes = new Geom_Plane(aPln); + break; + } + case GeomAbs_Cylinder: + { + gp_Cylinder aCyl; + if (aSampleSurf.GetType() == GeomAbs_Cylinder) + aCyl = aSampleSurf.Cylinder(); + if (aCanonRec.IsCylinder(tol, aCyl)) + aRes = new Geom_CylindricalSurface(aCyl); + break; + } + case GeomAbs_Cone: + { + gp_Cone aCon; + if (aSampleSurf.GetType() == GeomAbs_Cone) + aCon = aSampleSurf.Cone(); + if (aCanonRec.IsCone(tol, aCon)) + aRes = new Geom_ConicalSurface(aCon); + break; + } + case GeomAbs_Sphere: + { + gp_Sphere aSph; + if (aSampleSurf.GetType() == GeomAbs_Sphere) + aSph = aSampleSurf.Sphere(); + if (aCanonRec.IsSphere(tol, aSph)) + aRes = new Geom_SphericalSurface(aSph); + break; + } + default: + break; + } + + if (!aRes.IsNull()) + { + DrawTrSurf::Set(a[1], aRes); + di << "Gap = " << aCanonRec.GetGap() << "\n"; + } + else + { + di << "Cannot get required surface" << "\n"; + } + return 0; +} +//======================================================================= +//function : getanacurve +//purpose : +//======================================================================= + +Standard_Integer getanacurve(Draw_Interpretor& di, + Standard_Integer n, const char** a) +{ + if (n < 3) { + di << "Usage: \n"; + di << "getanacurve res shape [target [tol]] \n"; + di << "target is reqired type of curve and can be: lin, cir, ell \n"; + return 1; + } + TopoDS_Shape sh = DBRep::Get(a[2]); + if (sh.IsNull()) return 1; + TopAbs_ShapeEnum aShType = sh.ShapeType(); + if (!(aShType == TopAbs_WIRE || aShType == TopAbs_EDGE)) + { + di << "Wrong shape type, shape can be wire or or edge \n"; + return 1; + } + + GeomAbs_CurveType aTargets[] = { GeomAbs_Line, GeomAbs_Circle, GeomAbs_Ellipse }; + Standard_Integer icurv = 0; + if (n > 3) + { + if (strcmp(a[3],"lin") == 0) + icurv = 0; + else if (strcmp(a[3], "cir") == 0) + icurv = 1; + else if (strcmp(a[3], "ell") == 0) + icurv = 2; + } + + Standard_Real tol = Precision::Confusion(); + if (n > 4) + tol = Draw::Atof(a[4]); + + ShapeAnalysis_CanonicalRecognition aCanonRec(sh); + Handle(Geom_Curve) aRes; + switch (aTargets[icurv]) + { + case GeomAbs_Line: + { + gp_Lin aLin; + if (aCanonRec.IsLine(tol, aLin)) + aRes = new Geom_Line(aLin); + break; + } + case GeomAbs_Circle: + { + gp_Circ aCirc; + if (aCanonRec.IsCircle(tol, aCirc)) + aRes = new Geom_Circle(aCirc); + break; + } + case GeomAbs_Ellipse: + { + gp_Elips anElips; + if (aCanonRec.IsEllipse(tol, anElips)) + aRes = new Geom_Ellipse(anElips); + break; + } + default: + break; + } + + if (!aRes.IsNull()) + { + DrawTrSurf::Set(a[1], aRes); + di << "Gap = " << aCanonRec.GetGap() << "\n"; + } + else + { + di << "Cannot get required curve" << "\n"; + } + return 0; + +} //======================================================================= //function : InitCommands @@ -1018,4 +1223,8 @@ static Standard_Integer checkedge(Draw_Interpretor& di, Standard_Integer argc, c theCommands.Add("getareacontour","wire ",__FILE__, getareacontour, groupold); theCommands.Add ("checkselfintersection","wire [face]", __FILE__,checkselfintersection,g); theCommands.Add ("checkedge","edge [face]", __FILE__,checkedge,g); + theCommands.Add("getanasurf", "getanasurf res shape [target [tol [sample]]] ", __FILE__, getanasurf, g); + theCommands.Add("getanacurve", "getanacurve res shape [target [tol]]", __FILE__, getanacurve, g); + + } diff --git a/src/ShapeAnalysis/FILES b/src/ShapeAnalysis/FILES index 50a8c769d3..d9510dcbe9 100755 --- a/src/ShapeAnalysis/FILES +++ b/src/ShapeAnalysis/FILES @@ -44,3 +44,11 @@ ShapeAnalysis_WireOrder.cxx ShapeAnalysis_WireOrder.hxx ShapeAnalysis_WireVertex.cxx ShapeAnalysis_WireVertex.hxx +ShapeAnalysis_CanonicalRecognition.cxx +ShapeAnalysis_CanonicalRecognition.hxx +ShapeAnalysis_FuncSphereLSDist.cxx +ShapeAnalysis_FuncSphereLSDist.hxx +ShapeAnalysis_FuncCylinderLSDist.cxx +ShapeAnalysis_FuncCylinderLSDist.hxx +ShapeAnalysis_FuncConeLSDist.cxx +ShapeAnalysis_FuncConeLSDist.hxx diff --git a/src/ShapeAnalysis/ShapeAnalysis_CanonicalRecognition.cxx b/src/ShapeAnalysis/ShapeAnalysis_CanonicalRecognition.cxx new file mode 100644 index 0000000000..3abf8c12b5 --- /dev/null +++ b/src/ShapeAnalysis/ShapeAnalysis_CanonicalRecognition.cxx @@ -0,0 +1,1415 @@ +// Created on: 2000-01-20 +// Created by: data exchange team +// Copyright (c) 2000-2014 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software{ you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +//Declaration of static functions +static Standard_Integer GetNbPars(const GeomAbs_CurveType theTarget); +static Standard_Integer GetNbPars(const GeomAbs_SurfaceType theTarget); +static Standard_Boolean SetConicParameters(const GeomAbs_CurveType theTarget, const Handle(Geom_Curve)& theConic, + gp_Ax2& thePos, TColStd_Array1OfReal& theParams); +static Standard_Boolean CompareConicParams(const GeomAbs_CurveType theTarget, const Standard_Real theTol, + const gp_Ax2& theRefPos, const TColStd_Array1OfReal& theRefParams, + const gp_Ax2& thePos, const TColStd_Array1OfReal& theParams); +static Standard_Boolean SetSurfParams(const GeomAbs_SurfaceType theTarget, const Handle(Geom_Surface)& theElemSurf, + gp_Ax3& thePos, TColStd_Array1OfReal& theParams); +static Standard_Boolean CompareSurfParams(const GeomAbs_SurfaceType theTarget, const Standard_Real theTol, + const gp_Ax3& theRefPos, const TColStd_Array1OfReal& theRefParams, + const gp_Ax3& thePos, const TColStd_Array1OfReal& theParams); +static Standard_Real DeviationSurfParams(const GeomAbs_SurfaceType theTarget, + const gp_Ax3& theRefPos, const TColStd_Array1OfReal& theRefParams, + const gp_Ax3& thePos, const TColStd_Array1OfReal& theParams); +static Standard_Boolean GetSamplePoints(const TopoDS_Wire& theWire, + const Standard_Real theTol, const Standard_Integer theMaxNbInt, + Handle(TColgp_HArray1OfXYZ)& thePoints); +static Standard_Real GetLSGap(const Handle(TColgp_HArray1OfXYZ)& thePoints, const GeomAbs_SurfaceType theTarget, + const gp_Ax3& thePos, const TColStd_Array1OfReal& theParams); +static void FillSolverData(const GeomAbs_SurfaceType theTarget, + const gp_Ax3& thePos, const TColStd_Array1OfReal& theParams, + math_Vector& theStartPoint, + math_Vector& theFBnd, math_Vector& theLBnd, const Standard_Real theRelDev = 0.2); +static void SetCanonicParameters(const GeomAbs_SurfaceType theTarget, + const math_Vector& theSol, gp_Ax3& thePos, TColStd_Array1OfReal& theParams); + + +//======================================================================= +//function : Constructor +//purpose : +//======================================================================= + +ShapeAnalysis_CanonicalRecognition::ShapeAnalysis_CanonicalRecognition() : + mySType(TopAbs_SHAPE), myGap(-1.), myStatus(-1) +{ +} + +ShapeAnalysis_CanonicalRecognition::ShapeAnalysis_CanonicalRecognition(const TopoDS_Shape& theShape) : + mySType(TopAbs_SHAPE), myGap(-1.), myStatus(-1) +{ + Init(theShape); +} + +void ShapeAnalysis_CanonicalRecognition::SetShape(const TopoDS_Shape& theShape) +{ + Init(theShape); +} + +//======================================================================= +//function : Init +//purpose : +//======================================================================= + +void ShapeAnalysis_CanonicalRecognition::Init(const TopoDS_Shape& theShape) +{ + TopAbs_ShapeEnum aT = theShape.ShapeType(); + switch (aT) + { + case TopAbs_SHELL: + case TopAbs_FACE: + case TopAbs_WIRE: + case TopAbs_EDGE: + myShape = theShape; + mySType = aT; + myStatus = 0; + break; + default : + myStatus = 1; + break; + } +} + +//======================================================================= +//function : IsElementarySurf +//purpose : +//======================================================================= + +Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsElementarySurf(const GeomAbs_SurfaceType theTarget, + const Standard_Real theTol, + gp_Ax3& thePos, TColStd_Array1OfReal& theParams) +{ + if (myStatus != 0) + return Standard_False; + // + if (mySType == TopAbs_FACE) + { + Handle(Geom_Surface) anElemSurf = GetSurface(TopoDS::Face(myShape), + theTol, GeomConvert_Target, theTarget, myGap, + myStatus); + if (anElemSurf.IsNull()) + return Standard_False; + // + Standard_Boolean isOK = SetSurfParams(theTarget, anElemSurf, thePos, theParams); + if (!isOK) + { + myStatus = 1; + return Standard_False; + } + return Standard_True; + } + else if (mySType == TopAbs_SHELL) + { + Handle(Geom_Surface) anElemSurf = GetSurface(TopoDS::Shell(myShape), theTol, + GeomConvert_Target, theTarget, myGap, myStatus); + if (anElemSurf.IsNull()) + { + return Standard_False; + } + Standard_Boolean isOK = SetSurfParams(theTarget, anElemSurf, thePos, theParams); + if (!isOK) + { + myStatus = 1; + return Standard_False; + } + return Standard_True; + } + else if (mySType == TopAbs_EDGE) + { + Handle(Geom_Surface) anElemSurf = GetSurface(TopoDS::Edge(myShape), + theTol, + GeomConvert_Target, theTarget, + thePos, theParams, + myGap, myStatus); + + Standard_Boolean isOK = SetSurfParams(theTarget, anElemSurf, thePos, theParams); + if (!isOK) + { + myStatus = 1; + return Standard_False; + } + return Standard_True; + } + else if (mySType == TopAbs_WIRE) + { + Handle(Geom_Surface) anElemSurf = GetSurface(TopoDS::Wire(myShape), + theTol, + GeomConvert_Target, theTarget, + thePos, theParams, + myGap, myStatus); + + Standard_Boolean isOK = SetSurfParams(theTarget, anElemSurf, thePos, theParams); + if (!isOK) + { + myStatus = 1; + return Standard_False; + } + return Standard_True; + } + return Standard_False; +} + +//======================================================================= +//function : IsPlane +//purpose : +//======================================================================= + +Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsPlane(const Standard_Real theTol, gp_Pln& thePln) +{ + gp_Ax3 aPos = thePln.Position(); + TColStd_Array1OfReal aParams(1, 1); + // + GeomAbs_SurfaceType aTarget = GeomAbs_Plane; + if (IsElementarySurf(aTarget, theTol, aPos, aParams)) + { + thePln.SetPosition(aPos); + return Standard_True; + } + //Try to build plane for wire or edge + if (mySType == TopAbs_EDGE || mySType == TopAbs_WIRE) + { + BRepLib_FindSurface aFndSurf(myShape, theTol, Standard_True, Standard_False); + if (aFndSurf.Found()) + { + Handle(Geom_Plane) aPlane = Handle(Geom_Plane)::DownCast(aFndSurf.Surface()); + thePln = aPlane->Pln(); + myGap = aFndSurf.ToleranceReached(); + myStatus = 0; + return Standard_True; + } + else + myStatus = 1; + } + return Standard_False; +} + +//======================================================================= +//function : IsCylinder +//purpose : +//======================================================================= + +Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsCylinder(const Standard_Real theTol, gp_Cylinder& theCyl) +{ + gp_Ax3 aPos = theCyl.Position(); + TColStd_Array1OfReal aParams(1, 1); + aParams(1) = theCyl.Radius(); + // + GeomAbs_SurfaceType aTarget = GeomAbs_Cylinder; + if (IsElementarySurf(aTarget, theTol, aPos, aParams)) + { + theCyl.SetPosition(aPos); + theCyl.SetRadius(aParams(1)); + return Standard_True; + } + + if (aParams(1) > Precision::Infinite()) + { + //Sample cylinder does not seem to be set, least square method is not applicable. + return Standard_False; + } + if (myShape.ShapeType() == TopAbs_EDGE || myShape.ShapeType() == TopAbs_WIRE) + { + //Try to build surface by least square method; + TopoDS_Wire aWire; + if (myShape.ShapeType() == TopAbs_EDGE) + { + BRep_Builder aBB; + aBB.MakeWire(aWire); + aBB.Add(aWire, myShape); + } + else + { + aWire = TopoDS::Wire(myShape); + } + Standard_Boolean isDone = GetSurfaceByLS(aWire, theTol, aTarget, aPos, aParams, myGap, myStatus); + if (isDone) + { + theCyl.SetPosition(aPos); + theCyl.SetRadius(aParams(1)); + return Standard_True; + } + } + return Standard_False; +} + +//======================================================================= +//function : IsCone +//purpose : +//======================================================================= + +Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsCone(const Standard_Real theTol, gp_Cone& theCone) +{ + gp_Ax3 aPos = theCone.Position(); + TColStd_Array1OfReal aParams(1, 2); + aParams(1) = theCone.SemiAngle(); + aParams(2) = theCone.RefRadius(); + // + GeomAbs_SurfaceType aTarget = GeomAbs_Cone; + if (IsElementarySurf(aTarget, theTol, aPos, aParams)) + { + theCone.SetPosition(aPos); + theCone.SetSemiAngle(aParams(1)); + theCone.SetRadius(aParams(2)); + return Standard_True; + } + + + if (aParams(2) > Precision::Infinite()) + { + //Sample cone does not seem to be set, least square method is not applicable. + return Standard_False; + } + if (myShape.ShapeType() == TopAbs_EDGE || myShape.ShapeType() == TopAbs_WIRE) + { + //Try to build surface by least square method; + TopoDS_Wire aWire; + if (myShape.ShapeType() == TopAbs_EDGE) + { + BRep_Builder aBB; + aBB.MakeWire(aWire); + aBB.Add(aWire, myShape); + } + else + { + aWire = TopoDS::Wire(myShape); + } + Standard_Boolean isDone = GetSurfaceByLS(aWire, theTol, aTarget, aPos, aParams, myGap, myStatus); + if (isDone) + { + theCone.SetPosition(aPos); + theCone.SetSemiAngle(aParams(1)); + theCone.SetRadius(aParams(2)); + return Standard_True; + } + } + return Standard_False; +} + +//======================================================================= +//function : IsSphere +//purpose : +//======================================================================= + +Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsSphere(const Standard_Real theTol, gp_Sphere& theSphere) +{ + gp_Ax3 aPos = theSphere.Position(); + TColStd_Array1OfReal aParams(1, 1); + aParams(1) = theSphere.Radius(); + // + GeomAbs_SurfaceType aTarget = GeomAbs_Sphere; + if (IsElementarySurf(aTarget, theTol, aPos, aParams)) + { + theSphere.SetPosition(aPos); + theSphere.SetRadius(aParams(1)); + return Standard_True; + } + // + if (aParams(1) > Precision::Infinite()) + { + //Sample sphere does not seem to be set, least square method is not applicable. + return Standard_False; + } + if (myShape.ShapeType() == TopAbs_EDGE || myShape.ShapeType() == TopAbs_WIRE) + { + //Try to build surface by least square method; + TopoDS_Wire aWire; + if (myShape.ShapeType() == TopAbs_EDGE) + { + BRep_Builder aBB; + aBB.MakeWire(aWire); + aBB.Add(aWire, myShape); + } + else + { + aWire = TopoDS::Wire(myShape); + } + Standard_Boolean isDone = GetSurfaceByLS(aWire, theTol, aTarget, aPos, aParams, myGap, myStatus); + if (isDone) + { + theSphere.SetPosition(aPos); + theSphere.SetRadius(aParams(1)); + return Standard_True; + } + } + return Standard_False; +} + +//======================================================================= +//function : IsConic +//purpose : +//======================================================================= + +Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsConic(const GeomAbs_CurveType theTarget, + const Standard_Real theTol, + gp_Ax2& thePos, TColStd_Array1OfReal& theParams) +{ + if (myStatus != 0) + return Standard_False; + + if (mySType == TopAbs_EDGE) + { + Handle(Geom_Curve) aConic = GetCurve(TopoDS::Edge(myShape), theTol, + GeomConvert_Target, theTarget, myGap, myStatus); + + if (aConic.IsNull()) + return Standard_False; + + Standard_Boolean isOK = SetConicParameters(theTarget, aConic, thePos, theParams); + + if(!isOK) + { + myStatus = 1; + return Standard_False; + } + return Standard_True; + } + else if (mySType == TopAbs_WIRE) + { + TopoDS_Iterator anIter(myShape); + if (!anIter.More()) + { + myStatus = 1; + return Standard_False; + } + gp_Ax2 aPos; + TColStd_Array1OfReal aParams(1, theParams.Length()); + const TopoDS_Shape& anEdge = anIter.Value(); + + Handle(Geom_Curve) aConic = GetCurve(TopoDS::Edge(anEdge), theTol, + GeomConvert_Target, theTarget, myGap, myStatus); + if (aConic.IsNull()) + { + return Standard_False; + } + Standard_Boolean isOK = SetConicParameters(theTarget, aConic, thePos, theParams); + if (!isOK) + { + myStatus = 1; + return Standard_False; + } + if (!anIter.More()) + { + return Standard_True; + } + else + { + anIter.Next(); + } + for (; anIter.More(); anIter.Next()) + { + aConic = GetCurve(TopoDS::Edge(anIter.Value()), theTol, + GeomConvert_Target, theTarget, myGap, myStatus); + if (aConic.IsNull()) + { + return Standard_False; + } + isOK = SetConicParameters(theTarget, aConic, aPos, aParams); + isOK = CompareConicParams(theTarget, theTol, thePos, theParams, aPos, aParams); + + if (!isOK) + { + return Standard_False; + } + } + return Standard_True; + } + myStatus = 1; + return Standard_False; +} + + +//======================================================================= +//function : IsLine +//purpose : +//======================================================================= + +Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsLine(const Standard_Real theTol, gp_Lin& theLin) +{ + gp_Ax2 aPos; + TColStd_Array1OfReal aParams(1, 1); + + GeomAbs_CurveType aTarget = GeomAbs_Line; + Standard_Boolean isOK = IsConic(aTarget, theTol, aPos, aParams); + if (isOK) + { + theLin.SetPosition(aPos.Axis()); + } + return isOK; +} + +//======================================================================= +//function : IsCircle +//purpose : +//======================================================================= + +Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsCircle(const Standard_Real theTol, gp_Circ& theCirc) +{ + gp_Ax2 aPos; + TColStd_Array1OfReal aParams(1, 1); + + GeomAbs_CurveType aTarget = GeomAbs_Circle; + Standard_Boolean isOK = IsConic(aTarget, theTol, aPos, aParams); + if (isOK) + { + theCirc.SetPosition(aPos); + theCirc.SetRadius(aParams(1)); + } + return isOK; +} + +//======================================================================= +//function : IsEllipse +//purpose : +//======================================================================= + +Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsEllipse(const Standard_Real theTol, gp_Elips& theElips) +{ + gp_Ax2 aPos; + TColStd_Array1OfReal aParams(1, 2); + + GeomAbs_CurveType aTarget = GeomAbs_Ellipse; + Standard_Boolean isOK = IsConic(aTarget, theTol, aPos, aParams); + if (isOK) + { + theElips.SetPosition(aPos); + theElips.SetMajorRadius(aParams(1)); + theElips.SetMinorRadius(aParams(2)); + } + return isOK; +} + + +//======================================================================= +//function : GetSurface +//purpose : +//======================================================================= + +Handle(Geom_Surface) ShapeAnalysis_CanonicalRecognition::GetSurface(const TopoDS_Face& theFace, + const Standard_Real theTol, + const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget, + Standard_Real& theGap, Standard_Integer& theStatus) +{ + theStatus = 0; + TopLoc_Location aLoc; + const Handle(Geom_Surface)& aSurf = BRep_Tool::Surface(theFace, aLoc); + if (aSurf.IsNull()) + { + theStatus = 1; + return aSurf; + } + GeomConvert_SurfToAnaSurf aConv(aSurf); + aConv.SetConvType(theType); + aConv.SetTarget(theTarget); + Handle(Geom_Surface) anAnaSurf = aConv.ConvertToAnalytical(theTol); + if (anAnaSurf.IsNull()) + return anAnaSurf; + // + if (!aLoc.IsIdentity()) + anAnaSurf->Transform(aLoc.Transformation()); + // + theGap = aConv.Gap(); + return anAnaSurf; +} + +//======================================================================= +//function : GetSurface +//purpose : +//======================================================================= + +Handle(Geom_Surface) ShapeAnalysis_CanonicalRecognition::GetSurface(const TopoDS_Shell& theShell, + const Standard_Real theTol, + const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget, + Standard_Real& theGap, Standard_Integer& theStatus) +{ + Handle(Geom_Surface) anElemSurf; + TopoDS_Iterator anIter(theShell); + if (!anIter.More()) + { + theStatus = 1; + return anElemSurf; + } + gp_Ax3 aPos1; + TColStd_Array1OfReal aParams1(1, Max(1, GetNbPars(theTarget))); + const TopoDS_Shape& aFace = anIter.Value(); + + anElemSurf = GetSurface(TopoDS::Face(aFace), theTol, + theType, theTarget, theGap, theStatus); + if (anElemSurf.IsNull()) + { + return anElemSurf; + } + SetSurfParams(theTarget, anElemSurf, aPos1, aParams1); + if (!anIter.More()) + { + return anElemSurf; + } + else + { + anIter.Next(); + } + gp_Ax3 aPos; + TColStd_Array1OfReal aParams(1, Max(1, GetNbPars(theTarget))); + Standard_Boolean isOK; + for (; anIter.More(); anIter.Next()) + { + Handle(Geom_Surface) anElemSurf1 = GetSurface(TopoDS::Face(anIter.Value()), theTol, + theType, theTarget, theGap, theStatus); + if (anElemSurf1.IsNull()) + { + return anElemSurf1; + } + SetSurfParams(theTarget, anElemSurf, aPos, aParams); + isOK = CompareSurfParams(theTarget, theTol, aPos1, aParams1, aPos, aParams); + + if (!isOK) + { + anElemSurf.Nullify(); + return anElemSurf; + } + } + return anElemSurf; +} + +//======================================================================= +//function : GetSurface +//purpose : +//======================================================================= + +Handle(Geom_Surface) ShapeAnalysis_CanonicalRecognition::GetSurface(const TopoDS_Edge& theEdge, + const Standard_Real theTol, + const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget, + gp_Ax3& thePos, TColStd_Array1OfReal& theParams, + Standard_Real& theGap, Standard_Integer& theStatus) +{ + //Get surface list + NCollection_Vector aSurfs; + NCollection_Vector aGaps; + Standard_Integer j = 0; + for (;;) { + j++; + Handle(Geom_Surface) aSurf; + TopLoc_Location aLoc; + Handle(Geom2d_Curve) aPCurve; + Standard_Real ff, ll; + BRep_Tool::CurveOnSurface(theEdge, aPCurve, aSurf, aLoc, ff, ll, j); + if (aSurf.IsNull()) { + break; + } + GeomConvert_SurfToAnaSurf aConv(aSurf); + aConv.SetConvType(theType); + aConv.SetTarget(theTarget); + Handle(Geom_Surface) anAnaSurf = aConv.ConvertToAnalytical(theTol); + if (anAnaSurf.IsNull()) + continue; + // + if (!aLoc.IsIdentity()) + anAnaSurf->Transform(aLoc.Transformation()); + // + aGaps.Append(aConv.Gap()); + aSurfs.Append(anAnaSurf); + } + + if (aSurfs.Size() == 0) + { + theStatus = 1; + return NULL; + } + + gp_Ax3 aPos; + Standard_Integer aNbPars = Max(1, GetNbPars(theTarget)); + TColStd_Array1OfReal aParams(1, aNbPars); + + Standard_Integer ifit = -1, i; + Standard_Real aMinDev = RealLast(); + if (aSurfs.Size() == 1) + { + ifit = 0; + } + else + { + for (i = 0; i < aSurfs.Size(); ++i) + { + SetSurfParams(theTarget, aSurfs(i), aPos, aParams); + Standard_Real aDev = DeviationSurfParams(theTarget, thePos, theParams, aPos, aParams); + if (aDev < aMinDev) + { + aMinDev = aDev; + ifit = i; + } + } + } + if (ifit >= 0) + { + SetSurfParams(theTarget, aSurfs(ifit), thePos, theParams); + theGap = aGaps(ifit); + return aSurfs(ifit); + } + else + { + theStatus = 1; + return NULL; + } +} + +//======================================================================= +//function : GetSurface +//purpose : +//======================================================================= + +Handle(Geom_Surface) ShapeAnalysis_CanonicalRecognition::GetSurface(const TopoDS_Wire& theWire, + const Standard_Real theTol, + const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget, + gp_Ax3& thePos, TColStd_Array1OfReal& theParams, + Standard_Real& theGap, Standard_Integer& theStatus) +{ + //Get surface list + NCollection_Vector aSurfs; + NCollection_Vector aGaps; + + TopoDS_Iterator anIter(theWire); + if (!anIter.More()) + { + // Empty wire + theStatus = 1; + return NULL; + } + //First edge + const TopoDS_Edge& anEdge1 = TopoDS::Edge(anIter.Value()); + gp_Ax3 aPos1 = thePos; + Standard_Integer aNbPars = GetNbPars(theTarget); + TColStd_Array1OfReal aParams1(1, Max(1, aNbPars)); + Standard_Real aGap1; + Standard_Integer i; + for (i = 1; i <= aNbPars; ++i) + { + aParams1(i) = theParams(i); + } + Handle(Geom_Surface) anElemSurf1 = GetSurface(anEdge1, theTol, + theType, theTarget, aPos1, aParams1, aGap1, theStatus); + if (theStatus != 0 || anElemSurf1.IsNull()) + { + return NULL; + } + anIter.Next(); + for(; anIter.More(); anIter.Next()) + { + gp_Ax3 aPos = aPos1; + aNbPars = GetNbPars(theTarget); + TColStd_Array1OfReal aParams(1, Max(1, aNbPars)); + for (i = 1; i <= aNbPars; ++i) + { + aParams(i) = aParams1(i); + } + Standard_Real aGap; + const TopoDS_Edge& anEdge = TopoDS::Edge(anIter.Value()); + Handle(Geom_Surface) anElemSurf = GetSurface(anEdge, theTol, + theType, theTarget, aPos, aParams, aGap, theStatus); + if (theStatus != 0 || anElemSurf.IsNull()) + { + return NULL; + } + Standard_Boolean isOK = CompareSurfParams(theTarget, theTol, + aPos1, aParams1, aPos, aParams); + if (!isOK) + { + return NULL; + } + } + + SetSurfParams(theTarget, anElemSurf1, thePos, theParams); + theGap = aGap1; + return anElemSurf1; + +} +//======================================================================= +//function : GetSurfaceByLS +//purpose : +//======================================================================= + +Standard_Boolean ShapeAnalysis_CanonicalRecognition::GetSurfaceByLS(const TopoDS_Wire& theWire, + const Standard_Real theTol, + const GeomAbs_SurfaceType theTarget, + gp_Ax3& thePos, TColStd_Array1OfReal& theParams, + Standard_Real& theGap, Standard_Integer& theStatus) +{ + Handle(TColgp_HArray1OfXYZ) aPoints; + Standard_Integer aNbMaxInt = 100; + if (!GetSamplePoints(theWire, theTol, aNbMaxInt, aPoints)) + return Standard_False; + + theGap = GetLSGap(aPoints, theTarget, thePos, theParams); + if (theGap <= theTol) + { + theStatus = 0; + return Standard_True; + } + + Standard_Integer aNbVar = 0; + if (theTarget == GeomAbs_Sphere) + aNbVar = 4; + else if (theTarget == GeomAbs_Cylinder) + aNbVar = 4; + else if (theTarget == GeomAbs_Cone) + aNbVar = 5; + else + return Standard_False; + + math_Vector aFBnd(1, aNbVar), aLBnd(1, aNbVar), aStartPoint(1, aNbVar); + + Standard_Real aRelDev = 0.2; //Customer can set parameters of sample surface + // with relative precision about aRelDev. + // For example, if radius of sample surface is R, + // it means, that "exact" vaue is in interav + //[R - aRelDev*R, R + aRelDev*R]. This intrrval is set + // for R as boundary values for dptimization algo. + FillSolverData(theTarget, thePos, theParams, + aStartPoint, aFBnd, aLBnd, aRelDev); + + // + Standard_Real aTol = Precision::Confusion(); + math_MultipleVarFunction* aPFunc; + ShapeAnalysis_FuncSphereLSDist aFuncSph(aPoints); + ShapeAnalysis_FuncCylinderLSDist aFuncCyl(aPoints, thePos.Direction()); + ShapeAnalysis_FuncConeLSDist aFuncCon(aPoints, thePos.Direction()); + if (theTarget == GeomAbs_Sphere) + { + aPFunc = (math_MultipleVarFunction*)&aFuncSph; + } + else if (theTarget == GeomAbs_Cylinder) + { + aPFunc = (math_MultipleVarFunction*)&aFuncCyl; + } + else if (theTarget == GeomAbs_Cone) + { + aPFunc = (math_MultipleVarFunction*)&aFuncCon; + } + else + aPFunc = NULL; + // + math_Vector aSteps(1, aNbVar); + Standard_Integer aNbInt = 10; + Standard_Integer i; + for (i = 1; i <= aNbVar; ++i) + { + aSteps(i) = (aLBnd(i) - aFBnd(i)) / aNbInt; + } + math_PSO aGlobSolver(aPFunc, aFBnd, aLBnd, aSteps); + Standard_Real aLSDist; + aGlobSolver.Perform(aSteps, aLSDist, aStartPoint); + SetCanonicParameters(theTarget, aStartPoint, thePos, theParams); + + theGap = GetLSGap(aPoints, theTarget, thePos, theParams); + if (theGap <= theTol) + { + theStatus = 0; + return Standard_True; + } + // + math_Matrix aDirMatrix(1, aNbVar, 1, aNbVar, 0.0); + for (i = 1; i <= aNbVar; i++) + aDirMatrix(i, i) = 1.0; + + if (theTarget == GeomAbs_Cylinder || theTarget == GeomAbs_Cone) + { + //Set search direction for location to be perpendicular to axis to avoid + //seaching along axis + const gp_Dir aDir = thePos.Direction(); + gp_Pln aPln(thePos.Location(), aDir); + gp_Dir aUDir = aPln.Position().XDirection(); + gp_Dir aVDir = aPln.Position().YDirection(); + for (i = 1; i <= 3; ++i) + { + aDirMatrix(i, 1) = aUDir.Coord(i); + aDirMatrix(i, 2) = aVDir.Coord(i); + gp_Dir aUVDir(aUDir.XYZ() + aVDir.XYZ()); + aDirMatrix(i, 3) = aUVDir.Coord(i); + } + } + + math_Powell aSolver(*aPFunc, aTol); + aSolver.Perform(*aPFunc, aStartPoint, aDirMatrix); + + if (aSolver.IsDone()) + { + aSolver.Location(aStartPoint); + theStatus = 0; + SetCanonicParameters(theTarget, aStartPoint, thePos, theParams); + theGap = GetLSGap(aPoints, theTarget, thePos, theParams); + theStatus = 0; + if(theGap <= theTol) + return Standard_True; + } + else + theStatus = 1; + + return Standard_False; + +} + +//======================================================================= +//function : GetCurve +//purpose : +//======================================================================= + +Handle(Geom_Curve) ShapeAnalysis_CanonicalRecognition::GetCurve(const TopoDS_Edge& theEdge, + const Standard_Real theTol, + const GeomConvert_ConvType theType, const GeomAbs_CurveType theTarget, + Standard_Real& theGap, Standard_Integer& theStatus) +{ + theStatus = 0; + TopLoc_Location aLoc; + Standard_Real f, l, nf, nl; + const Handle(Geom_Curve)& aCurv = BRep_Tool::Curve(theEdge, aLoc, f, l); + if (aCurv.IsNull()) + { + theStatus = 1; + return aCurv; + } + GeomConvert_CurveToAnaCurve aConv(aCurv); + aConv.SetConvType(theType); + aConv.SetTarget(theTarget); + Handle(Geom_Curve) anAnaCurv; + aConv.ConvertToAnalytical(theTol, anAnaCurv, f, l, nf, nl); + if (anAnaCurv.IsNull()) + return anAnaCurv; + // + if (!aLoc.IsIdentity()) + anAnaCurv->Transform(aLoc.Transformation()); + // + theGap = aConv.Gap(); + return anAnaCurv; +} + +//Static methods +//======================================================================= +//function : GetNbPars +//purpose : +//======================================================================= + +Standard_Integer GetNbPars(const GeomAbs_CurveType theTarget) +{ + Standard_Integer aNbPars = 0; + switch (theTarget) + { + case GeomAbs_Line: + aNbPars = 0; + break; + case GeomAbs_Circle: + aNbPars = 1; + break; + case GeomAbs_Ellipse: + aNbPars = 2; + break; + default: + aNbPars = 0; + break; + } + + return aNbPars; +} + +//======================================================================= +//function : GetNbPars +//purpose : +//======================================================================= + +Standard_Integer GetNbPars(const GeomAbs_SurfaceType theTarget) +{ + Standard_Integer aNbPars = 0; + switch (theTarget) + { + case GeomAbs_Plane: + aNbPars = 0; + break; + case GeomAbs_Cylinder: + case GeomAbs_Sphere: + aNbPars = 1; + break; + case GeomAbs_Cone: + aNbPars = 2; + break; + default: + aNbPars = 0; + break; + } + + return aNbPars; +} + +//======================================================================= +//function : SetConicParameters +//purpose : +//======================================================================= + +Standard_Boolean SetConicParameters(const GeomAbs_CurveType theTarget, const Handle(Geom_Curve)& theConic, + gp_Ax2& thePos, TColStd_Array1OfReal& theParams) +{ + if (theConic.IsNull()) + return Standard_False; + GeomAdaptor_Curve aGAC(theConic); + if(aGAC.GetType() != theTarget) + return Standard_False; + + if (theTarget == GeomAbs_Line) + { + gp_Lin aLin = aGAC.Line(); + thePos.SetAxis(aLin.Position()); + } + else if (theTarget == GeomAbs_Circle) + { + gp_Circ aCirc = aGAC.Circle(); + thePos = aCirc.Position(); + theParams(1) = aCirc.Radius(); + } + else if (theTarget == GeomAbs_Ellipse) + { + gp_Elips anElips = aGAC.Ellipse(); + thePos = anElips.Position(); + theParams(1) = anElips.MajorRadius(); + theParams(2) = anElips.MinorRadius(); + } + else + return Standard_False; + return Standard_True; + +} + +//======================================================================= +//function : CompareConicParams +//purpose : +//======================================================================= + +Standard_Boolean CompareConicParams(const GeomAbs_CurveType theTarget, const Standard_Real theTol, + const gp_Ax2& theRefPos, const TColStd_Array1OfReal& theRefParams, + const gp_Ax2& thePos, const TColStd_Array1OfReal& theParams) +{ + Standard_Integer i, aNbPars = GetNbPars(theTarget); + + for (i = 1; i <= aNbPars; ++i) + { + if (Abs(theRefParams(i) - theParams(i)) > theTol) + return Standard_False; + } + + Standard_Real anAngTol = theTol / (2. * M_PI); + Standard_Real aTol = theTol; + if (theTarget == GeomAbs_Line) + aTol = Precision::Infinite(); + + const gp_Ax1& aRef = theRefPos.Axis(); + const gp_Ax1& anAx1 = thePos.Axis(); + gp_Ax1 anAx1Rev = anAx1.Reversed(); + + if (aRef.IsCoaxial(anAx1, anAngTol, aTol) || aRef.IsCoaxial(anAx1Rev, anAngTol, aTol)) + { + return Standard_True; + } + return Standard_False; +} + +//======================================================================= +//function : SetSurfParams +//purpose : +//======================================================================= + +Standard_Boolean SetSurfParams(const GeomAbs_SurfaceType theTarget, const Handle(Geom_Surface)& theElemSurf, + gp_Ax3& thePos, TColStd_Array1OfReal& theParams) +{ + // + if (theElemSurf.IsNull()) + return Standard_False; + GeomAdaptor_Surface aGAS(theElemSurf); + if (aGAS.GetType() != theTarget) + return Standard_False; + + Standard_Real aNbPars = GetNbPars(theTarget); + if (theParams.Length() < aNbPars) + return Standard_False; + + if (theTarget == GeomAbs_Plane) + { + gp_Pln aPln = aGAS.Plane(); + thePos = aPln.Position(); + } + else if (theTarget == GeomAbs_Cylinder) + { + gp_Cylinder aCyl = aGAS.Cylinder(); + thePos = aCyl.Position(); + theParams(1) = aCyl.Radius(); + } + else if (theTarget == GeomAbs_Cone) + { + gp_Cone aCon = aGAS.Cone(); + thePos = aCon.Position(); + theParams(1) = aCon.SemiAngle(); + theParams(2) = aCon.RefRadius(); + } + else if (theTarget == GeomAbs_Sphere) + { + gp_Sphere aSph = aGAS.Sphere(); + thePos = aSph.Position(); + theParams(1) = aSph.Radius(); + } + else + { + return Standard_False; + } + return Standard_True; +} + +//======================================================================= +//function : CompareSurfParams +//purpose : +//======================================================================= + +Standard_Boolean CompareSurfParams(const GeomAbs_SurfaceType theTarget, const Standard_Real theTol, + const gp_Ax3& theRefPos, const TColStd_Array1OfReal& theRefParams, + const gp_Ax3& thePos, const TColStd_Array1OfReal& theParams) +{ + if (theTarget != GeomAbs_Plane) + { + if (Abs(theRefParams(1) - theParams(1)) > theTol) + { + return Standard_False; + } + } + // + if (theTarget == GeomAbs_Sphere) + { + gp_Pnt aRefLoc = theRefPos.Location(), aLoc = thePos.Location(); + if (aRefLoc.SquareDistance(aLoc) <= theTol*theTol) + { + return Standard_True; + } + return Standard_False; + } + // + Standard_Real anAngTol = theTol / (2. * M_PI); + Standard_Real aTol = theTol; + if (theTarget == GeomAbs_Cylinder || theTarget == GeomAbs_Cone) + { + aTol = Precision::Infinite(); + } + + const gp_Ax1& aRef = theRefPos.Axis(); + const gp_Ax1& anAx1 = thePos.Axis(); + gp_Ax1 anAx1Rev = anAx1.Reversed(); + if (!(aRef.IsCoaxial(anAx1, anAngTol, aTol) || aRef.IsCoaxial(anAx1Rev, anAngTol, aTol))) + { + return Standard_False; + } + + if (theTarget == GeomAbs_Cone) + { + gp_Cone aRefCone(theRefPos, theRefParams(1), theRefParams(2)); + gp_Cone aCone(thePos, theParams(1), theParams(2)); + gp_Pnt aRefApex = aRefCone.Apex(); + gp_Pnt anApex = aCone.Apex(); + if (aRefApex.SquareDistance(anApex) <= theTol*theTol) + { + return Standard_True; + } + return Standard_False; + } + + return Standard_True; +} + +//======================================================================= +//function : DeviationSurfParams +//purpose : +//======================================================================= + +Standard_Real DeviationSurfParams(const GeomAbs_SurfaceType theTarget, + const gp_Ax3& theRefPos, const TColStd_Array1OfReal& theRefParams, + const gp_Ax3& thePos, const TColStd_Array1OfReal& theParams) +{ + Standard_Real aDevPars = 0.; + if (theTarget != GeomAbs_Plane) + { + aDevPars = Abs(theRefParams(1) - theParams(1)); + } + // + if (theTarget == GeomAbs_Sphere) + { + gp_Pnt aRefLoc = theRefPos.Location(), aLoc = thePos.Location(); + aDevPars += aRefLoc.Distance(aLoc); + } + else + { + const gp_Dir& aRefDir = theRefPos.Direction(); + const gp_Dir& aDir = thePos.Direction(); + Standard_Real anAngDev = (1. - Abs(aRefDir * aDir)); + aDevPars += anAngDev; + } + + return aDevPars; +} + +//======================================================================= +//function : GetSamplePoints +//purpose : +//======================================================================= + +Standard_Boolean GetSamplePoints(const TopoDS_Wire& theWire, + const Standard_Real theTol, + const Standard_Integer theMaxNbInt, + Handle(TColgp_HArray1OfXYZ)& thePoints) +{ + NCollection_Vector aLengths; + NCollection_Vector aCurves; + NCollection_Vector aPoints; + Standard_Real aTol = Max(1.e-3, theTol/10.); + Standard_Real aTotalLength = 0.; + TopoDS_Iterator anEIter(theWire); + for (; anEIter.More(); anEIter.Next()) + { + const TopoDS_Edge& anE = TopoDS::Edge(anEIter.Value()); + if (BRep_Tool::Degenerated(anE)) + continue; + BRepAdaptor_Curve aBAC(anE); + Standard_Real aClength = GCPnts_AbscissaPoint::Length(aBAC, aTol); + aTotalLength += aClength; + aCurves.Append(aBAC); + aLengths.Append(aClength); + } + + if(aTotalLength < theTol) + return Standard_False; + + Standard_Integer i, aNb = aLengths.Length(); + for (i = 0; i < aNb; ++i) + { + const BRepAdaptor_Curve& aC = aCurves(i); + Standard_Real aClength = GCPnts_AbscissaPoint::Length(aC, aTol); + Standard_Integer aNbPoints = RealToInt(aClength / aTotalLength * theMaxNbInt + 1); + aNbPoints = Max(2, aNbPoints); + GCPnts_QuasiUniformAbscissa aPointGen(aC, aNbPoints); + if (!aPointGen.IsDone()) + continue; + aNbPoints = aPointGen.NbPoints(); + Standard_Integer j; + for (j = 1; j <= aNbPoints; ++j) + { + Standard_Real t = aPointGen.Parameter(j); + gp_Pnt aP = aC.Value(t); + aPoints.Append(aP.XYZ()); + } + } + + if (aPoints.Length() < 1) + return Standard_False; + + thePoints = new TColgp_HArray1OfXYZ(1, aPoints.Length()); + for (i = 0; i < aPoints.Length(); ++i) + { + thePoints->SetValue(i + 1, aPoints(i)); + } + + return Standard_True; + +} + +//======================================================================= +//function : GetLSGap +//purpose : +//======================================================================= + +static Standard_Real GetLSGap(const Handle(TColgp_HArray1OfXYZ)& thePoints, const GeomAbs_SurfaceType theTarget, + const gp_Ax3& thePos, const TColStd_Array1OfReal& theParams) +{ + + Standard_Real aGap = 0.; + gp_XYZ aLoc = thePos.Location().XYZ(); + gp_Vec aDir(thePos.Direction()); + + + Standard_Integer i; + if (theTarget == GeomAbs_Sphere) + { + Standard_Real anR = theParams(1); + for (i = thePoints->Lower(); i <= thePoints->Upper(); ++i) + { + gp_XYZ aD = thePoints->Value(i) - aLoc; + aGap = Max(aGap, Abs((aD.Modulus() - anR))); + } + } + else if (theTarget == GeomAbs_Cylinder) + { + Standard_Real anR = theParams(1); + for (i = thePoints->Lower(); i <= thePoints->Upper(); ++i) + { + gp_Vec aD(thePoints->Value(i) - aLoc); + aD.Cross(aDir); + aGap = Max(aGap, Abs((aD.Magnitude() - anR))); + } + } + else if (theTarget == GeomAbs_Cone) + { + Standard_Real anAng = theParams(1); + Standard_Real anR = theParams(2); + for (i = thePoints->Lower(); i <= thePoints->Upper(); ++i) + { + Standard_Real u, v; + gp_Pnt aPi(thePoints->Value(i)); + ElSLib::ConeParameters(thePos, anR, anAng, aPi, u, v); + gp_Pnt aPp; + ElSLib::ConeD0(u, v, thePos, anR, anAng, aPp); + aGap = Max(aGap, aPi.SquareDistance(aPp)); + } + aGap = Sqrt(aGap); + } + + return aGap; +} + +//======================================================================= +//function : FillSolverData +//purpose : +//======================================================================= + +void FillSolverData(const GeomAbs_SurfaceType theTarget, + const gp_Ax3& thePos, const TColStd_Array1OfReal& theParams, + math_Vector& theStartPoint, + math_Vector& theFBnd, math_Vector& theLBnd, const Standard_Real theRelDev) +{ + if (theTarget == GeomAbs_Sphere || theTarget == GeomAbs_Cylinder) + { + theStartPoint(1) = thePos.Location().X(); + theStartPoint(2) = thePos.Location().Y(); + theStartPoint(3) = thePos.Location().Z(); + theStartPoint(4) = theParams(1); + Standard_Real aDR = theRelDev * theParams(1); + Standard_Real aDXYZ = aDR; + Standard_Integer i; + for (i = 1; i <= 3; ++i) + { + theFBnd(i) = theStartPoint(i) - aDXYZ; + theLBnd(i) = theStartPoint(i) + aDXYZ; + } + theFBnd(4) = theStartPoint(4) - aDR; + theLBnd(4) = theStartPoint(4) + aDR; + } + if (theTarget == GeomAbs_Cone) + { + theStartPoint(1) = thePos.Location().X(); + theStartPoint(2) = thePos.Location().Y(); + theStartPoint(3) = thePos.Location().Z(); + theStartPoint(4) = theParams(1); //SemiAngle + theStartPoint(5) = theParams(2); //Radius + Standard_Real aDR = theRelDev * theParams(2); + if (aDR < Precision::Confusion()) + { + aDR = 0.1; + } + Standard_Real aDXYZ = aDR; + Standard_Real aDAng = theRelDev * Abs(theParams(1)); + Standard_Integer i; + for (i = 1; i <= 3; ++i) + { + theFBnd(i) = theStartPoint(i) - aDXYZ; + theLBnd(i) = theStartPoint(i) + aDXYZ; + } + if (theParams(1) >= 0.) + { + theFBnd(4) = theStartPoint(4) - aDAng; + theLBnd(4) = Min(M_PI_2, theStartPoint(4) + aDR); + } + else + { + theFBnd(4) = Max(-M_PI_2, theStartPoint(4) - aDAng); + theLBnd(4) = theStartPoint(4) + aDAng; + } + theFBnd(5) = theStartPoint(5) - aDR; + theLBnd(5) = theStartPoint(5) + aDR; + } + +} + +//======================================================================= +//function : SetCanonicParameters +//purpose : +//======================================================================= + +void SetCanonicParameters(const GeomAbs_SurfaceType theTarget, + const math_Vector& theSol, gp_Ax3& thePos, TColStd_Array1OfReal& theParams) +{ + gp_Pnt aLoc(theSol(1), theSol(2), theSol(3)); + thePos.SetLocation(aLoc); + if (theTarget == GeomAbs_Sphere || theTarget == GeomAbs_Cylinder) + { + theParams(1) = theSol(4);//radius + } + else if (theTarget == GeomAbs_Cone) + { + theParams(1) = theSol(4);//semiangle + theParams(2) = theSol(5);//radius + } + +} \ No newline at end of file diff --git a/src/ShapeAnalysis/ShapeAnalysis_CanonicalRecognition.hxx b/src/ShapeAnalysis/ShapeAnalysis_CanonicalRecognition.hxx new file mode 100644 index 0000000000..7fa079b82c --- /dev/null +++ b/src/ShapeAnalysis/ShapeAnalysis_CanonicalRecognition.hxx @@ -0,0 +1,169 @@ +// Copyright (c) 2022 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _ShapeAnalysis_CanonicalRecognition_HeaderFile +#define _ShapeAnalysis_CanonicalRecognition_HeaderFile + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +class gp_Pln; +class gp_Cone; +class gp_Cylinder; +class gp_Sphere; +class gp_Lin; +class gp_Circ; +class gp_Elips; +class Geom_Curve; +class Geom_Surface; + +//! This class provides operators for analysis surfaces and curves of shapes +//! in order to find out more simple geometry entities, which could replace +//! existing complex (for exampe, BSpline) geometry objects with given tolerance. +class ShapeAnalysis_CanonicalRecognition +{ +public: + + DEFINE_STANDARD_ALLOC + + + //! Empty constructor + Standard_EXPORT ShapeAnalysis_CanonicalRecognition(); + + //! constructor with shape initialisation + Standard_EXPORT ShapeAnalysis_CanonicalRecognition(const TopoDS_Shape& theShape); + + //! Sets shape + Standard_EXPORT void SetShape(const TopoDS_Shape& theShape); + + //! Returns input shape + const TopoDS_Shape& GetShape() const + { + return myShape; + } + + //! Returns deviation between input geometry entity and analytical entity + Standard_Real GetGap() const + { + return myGap; + } + + //! Returns status of operation. + //! Current meaning of possible values of status: + //! -1 - algorithm is not initalazed by shape + //! 0 - no errors + //! 1 - error during any operation (usually - because of wrong input data) + //! Any operation (calling any methods like IsPlane(...), ...) can be performed + //! when current staue is equal 0. + //! If after any operation status != 0, it is necessary to set it 0 by method ClearStatus() + //! before calling other operation. + Standard_Integer GetStatus() const + { + return myStatus; + } + + //! Returns status to be equal 0. + void ClearStatus() + { + myStatus = 0; + } + + //! Returns true if the underlined surface can be represent by plane with tolerance theTol + //! and sets in thePln the result plane. + Standard_EXPORT Standard_Boolean IsPlane(const Standard_Real theTol, gp_Pln& thePln); + + //! Returns true if the underlined surface can be represent by cylindrical one with tolerance theTol + //! and sets in theCyl the result cylinrical surface. + Standard_EXPORT Standard_Boolean IsCylinder(const Standard_Real theTol, gp_Cylinder& theCyl); + + //! Returns true if the underlined surface can be represent by conical one with tolerance theTol + //! and sets in theCone the result conical surface. + Standard_EXPORT Standard_Boolean IsCone(const Standard_Real theTol, gp_Cone& theCone); + + //! Returns true if the underlined surface can be represent by spherical one with tolerance theTol + //! and sets in theSphere the result spherical surface. + Standard_EXPORT Standard_Boolean IsSphere(const Standard_Real theTol, gp_Sphere& theSphere); + + //! Returns true if the underlined curve can be represent by line with tolerance theTol + //! and sets in theLin the result line. + Standard_EXPORT Standard_Boolean IsLine(const Standard_Real theTol, gp_Lin& theLin); + + //! Returns true if the underlined curve can be represent by circle with tolerance theTol + //! and sets in theCirc the result circle. + Standard_EXPORT Standard_Boolean IsCircle(const Standard_Real theTol, gp_Circ& theCirc); + + //! Returns true if the underlined curve can be represent by ellipse with tolerance theTol + //! and sets in theCirc the result ellipse. + Standard_EXPORT Standard_Boolean IsEllipse(const Standard_Real theTol, gp_Elips& theElips); + + +private: + Standard_Boolean IsElementarySurf(const GeomAbs_SurfaceType theTarget, const Standard_Real theTol, + gp_Ax3& thePos, TColStd_Array1OfReal& theParams); + + Standard_Boolean IsConic(const GeomAbs_CurveType theTarget, const Standard_Real theTol, + gp_Ax2& thePos, TColStd_Array1OfReal& theParams); + + static Handle(Geom_Surface) GetSurface(const TopoDS_Face& theFace, const Standard_Real theTol, + const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget, + Standard_Real& theGap, Standard_Integer& theStatus); + + static Handle(Geom_Surface) GetSurface(const TopoDS_Shell& theShell, const Standard_Real theTol, + const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget, + Standard_Real& theGap, Standard_Integer& theStatus); + + static Handle(Geom_Surface) GetSurface(const TopoDS_Edge& theEdge, const Standard_Real theTol, + const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget, + gp_Ax3& thePos, TColStd_Array1OfReal& theParams, + Standard_Real& theGap, Standard_Integer& theStatus); + + static Handle(Geom_Surface) GetSurface(const TopoDS_Wire& theWire, const Standard_Real theTol, + const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget, + gp_Ax3& thePos, TColStd_Array1OfReal& theParams, + Standard_Real& theGap, Standard_Integer& theStatus); + + static Handle(Geom_Curve) GetCurve(const TopoDS_Edge& theEdge, const Standard_Real theTol, + const GeomConvert_ConvType theType, const GeomAbs_CurveType theTarget, + Standard_Real& theGap, Standard_Integer& theStatus); + + static Standard_Boolean GetSurfaceByLS(const TopoDS_Wire& theWire, const Standard_Real theTol, + const GeomAbs_SurfaceType theTarget, + gp_Ax3& thePos, TColStd_Array1OfReal& theParams, + Standard_Real& theGap, Standard_Integer& theStatus); + + void Init(const TopoDS_Shape& theShape); + +private: + + TopoDS_Shape myShape; + TopAbs_ShapeEnum mySType; + Standard_Real myGap; + Standard_Integer myStatus; + +}; + +#endif // _ShapeAnalysis_CanonicalRecognition_HeaderFile diff --git a/src/ShapeAnalysis/ShapeAnalysis_FuncConeLSDist.cxx b/src/ShapeAnalysis/ShapeAnalysis_FuncConeLSDist.cxx new file mode 100644 index 0000000000..865d9e872a --- /dev/null +++ b/src/ShapeAnalysis/ShapeAnalysis_FuncConeLSDist.cxx @@ -0,0 +1,68 @@ +// Copyright (c) 1995-1999 Matra Datavision +// Copyright (c) 1999-2022 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + + +#include +#include +#include +#include +#include +#include + +//======================================================================= +//function : ShapeAnalysis_FuncConeLSDist +//purpose : +//======================================================================= +ShapeAnalysis_FuncConeLSDist::ShapeAnalysis_FuncConeLSDist( + const Handle(TColgp_HArray1OfXYZ)& thePoints, + const gp_Dir& theDir): + myPoints(thePoints), myDir(theDir) +{ +} + +//======================================================================= +//function : NbVariables +//purpose : +//======================================================================= +Standard_Integer ShapeAnalysis_FuncConeLSDist::NbVariables () const +{ + return 5; +} + +//======================================================================= +//function : Value +//purpose : +//======================================================================= +Standard_Boolean ShapeAnalysis_FuncConeLSDist::Value(const math_Vector& X, Standard_Real& F) +{ + gp_Pnt aLoc(X(1), X(2), X(3)); + Standard_Real aSemiAngle = X(4), anR = X(5); + gp_Ax3 aPos(aLoc, myDir); + + F = 0.; + Standard_Integer i; + for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i) + { + Standard_Real u, v; + gp_Pnt aPi(myPoints->Value(i)); + ElSLib::ConeParameters(aPos, anR, aSemiAngle, aPi, u, v); + gp_Pnt aPp; + ElSLib::ConeD0(u, v, aPos, anR, aSemiAngle, aPp); + F += aPi.SquareDistance(aPp); + } + + return Standard_True; +} + + diff --git a/src/ShapeAnalysis/ShapeAnalysis_FuncConeLSDist.hxx b/src/ShapeAnalysis/ShapeAnalysis_FuncConeLSDist.hxx new file mode 100644 index 0000000000..381b86cc3f --- /dev/null +++ b/src/ShapeAnalysis/ShapeAnalysis_FuncConeLSDist.hxx @@ -0,0 +1,66 @@ +// Copyright (c) 1991-1999 Matra Datavision +// Copyright (c) 1999-2022 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _ShapeAnalysis_FuncConeLSDist_HeaderFile +#define _ShapeAnalysis_FuncConeLSDist_HeaderFile + +#include +#include + +#include +#include +#include +#include + +//! Function for search of Cone canonic parameters: coordinates of center local coordinate system, +//! direction of axis, radius and semi-angle from set of points +//! by least square method. +//! +//! +class ShapeAnalysis_FuncConeLSDist : public math_MultipleVarFunction +{ +public: + + DEFINE_STANDARD_ALLOC + + //! Constructor. + Standard_EXPORT ShapeAnalysis_FuncConeLSDist() {}; + + Standard_EXPORT ShapeAnalysis_FuncConeLSDist(const Handle(TColgp_HArray1OfXYZ)& thePoints, + const gp_Dir& theDir); + + void SetPoints(const Handle(TColgp_HArray1OfXYZ)& thePoints) + { + myPoints = thePoints; + } + + void SetDir(const gp_Dir& theDir) + { + myDir = theDir; + } + + //! Number of variables. + Standard_EXPORT Standard_Integer NbVariables() const Standard_OVERRIDE; + + //! Value. + Standard_EXPORT Standard_Boolean Value(const math_Vector& X,Standard_Real& F) Standard_OVERRIDE; + + +private: + + Handle(TColgp_HArray1OfXYZ) myPoints; + gp_Dir myDir; + +}; +#endif // _ShapeAnalysis_FuncConeLSDist_HeaderFile diff --git a/src/ShapeAnalysis/ShapeAnalysis_FuncCylinderLSDist.cxx b/src/ShapeAnalysis/ShapeAnalysis_FuncCylinderLSDist.cxx new file mode 100644 index 0000000000..681bef12f2 --- /dev/null +++ b/src/ShapeAnalysis/ShapeAnalysis_FuncCylinderLSDist.cxx @@ -0,0 +1,140 @@ +// Copyright (c) 1995-1999 Matra Datavision +// Copyright (c) 1999-2022 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + + +#include +#include +#include +#include + +//======================================================================= +//function : ShapeAnalysis_FuncCylinderLSDist +//purpose : +//======================================================================= +ShapeAnalysis_FuncCylinderLSDist::ShapeAnalysis_FuncCylinderLSDist( + const Handle(TColgp_HArray1OfXYZ)& thePoints, + const gp_Dir& theDir): + myPoints(thePoints), myDir(theDir) +{ +} + +//======================================================================= +//function : NbVariables +//purpose : +//======================================================================= +Standard_Integer ShapeAnalysis_FuncCylinderLSDist::NbVariables () const +{ + return 4; +} + +//======================================================================= +//function : Value +//purpose : +//======================================================================= +Standard_Boolean ShapeAnalysis_FuncCylinderLSDist::Value(const math_Vector& X,Standard_Real& F) +{ + gp_XYZ aLoc(X(1), X(2), X(3)); + Standard_Real anR2 = X(4)*X(4); + + F = 0.; + Standard_Integer i; + for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i) + { + gp_Vec aV(myPoints->Value(i) - aLoc); + Standard_Real aD2 = aV.CrossSquareMagnitude(myDir); + Standard_Real d = aD2 - anR2; + F += d * d; + } + + return Standard_True; +} + +//======================================================================= +//function : Gradient +//purpose : +//======================================================================= +Standard_Boolean ShapeAnalysis_FuncCylinderLSDist::Gradient(const math_Vector& X,math_Vector& G) + +{ + gp_XYZ aLoc(X(1), X(2), X(3)); + Standard_Real anR = X(4), anR2 = anR * anR; + Standard_Real x = myDir.X(), y = myDir.Y(), z = myDir.Z(); + G.Init(0.); + + Standard_Integer i; + for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i) + { + gp_Vec aV(myPoints->Value(i) - aLoc); + Standard_Real aD2 = aV.CrossSquareMagnitude(myDir); + Standard_Real d = aD2 - anR2; + Standard_Real Dx0 = 2.*(aV.Z()*x - aV.X()*z)*z + -2.*(aV.X()*y - aV.Y()*x)*y; + Standard_Real Dy0 = -2.*(aV.Y()*z - aV.Z()*y)*z + +2.*(aV.X()*y - aV.Y()*x)*x; + Standard_Real Dz0 = 2.*(aV.Y()*z - aV.Z()*y)*y + -2.*(aV.Z()*x - aV.X()*z)*x; + + G(1) += d * Dx0; + G(2) += d * Dy0; + G(3) += d * Dz0; + // + G(4) += d; + } + + G *= 2; + G(6) *= -2.*anR; + + return Standard_True; +} + +//======================================================================= +//function : Values +//purpose : +//======================================================================= +Standard_Boolean ShapeAnalysis_FuncCylinderLSDist::Values(const math_Vector& X,Standard_Real& F,math_Vector& G) +{ + gp_XYZ aLoc(X(1), X(2), X(3)); + Standard_Real anR = X(4), anR2 = anR * anR; + Standard_Real x = myDir.X(), y = myDir.Y(), z = myDir.Z(); + + F = 0.; + G.Init(0.); + Standard_Integer i; + for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i) + { + gp_Vec aV(myPoints->Value(i) - aLoc); + Standard_Real aD2 = aV.CrossSquareMagnitude(myDir); + Standard_Real d = aD2 - anR2; + Standard_Real Dx0 = 2.*(aV.Z()*x - aV.X()*z)*z + - 2.*(aV.X()*y - aV.Y()*x)*y; + Standard_Real Dy0 = -2.*(aV.Y()*z - aV.Z()*y)*z + + 2.*(aV.X()*y - aV.Y()*x)*x; + Standard_Real Dz0 = 2.*(aV.Y()*z - aV.Z()*y)*y + - 2.*(aV.Z()*x - aV.X()*z)*x; + + G(1) += d * Dx0; + G(2) += d * Dy0; + G(3) += d * Dz0; + // + G(4) += d; + // + F += d * d; + } + + G *= 2; + G(4) *= -2.*anR; + + return true; +} + diff --git a/src/ShapeAnalysis/ShapeAnalysis_FuncCylinderLSDist.hxx b/src/ShapeAnalysis/ShapeAnalysis_FuncCylinderLSDist.hxx new file mode 100644 index 0000000000..3375e522e7 --- /dev/null +++ b/src/ShapeAnalysis/ShapeAnalysis_FuncCylinderLSDist.hxx @@ -0,0 +1,100 @@ +// Copyright (c) 1991-1999 Matra Datavision +// Copyright (c) 1999-2022 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _ShapeAnalysis_FuncCylinderLSDist_HeaderFile +#define _ShapeAnalysis_FuncCylinderLSDist_HeaderFile + +#include +#include + +#include +#include +#include +#include + +//! Function for search of cylinder canonic parameters: coordinates of center local coordinate system, +//! direction of axis and radius from set of points +//! by least square method. +//! +//! The class inherits math_MultipleVarFunctionWithGradient and thus is intended +//! for use in math_BFGS algorithm. +//! +//! Parametrisation: +//! Cylinder is defined by its axis and radius. Axis is defined by 3 cartesian coordinats it location x0, y0, z0 +//! and direction, which is constant and set by user: +//! dir.x, dir.y, dir.z +//! The criteria is: +//! F(x0, y0, z0, theta, phi, R) = Sum[|(P(i) - Loc)^dir|^2 - R^2]^2 => min +//! P(i) is i-th sample point, Loc, dir - axis location and direction, R - radius +//! +//! The square vector product |(P(i) - Loc)^dir|^2 is: +//! +//! [(y - y0)*dir.z - (z - z0)*dir.y]^2 + +//! [(z - z0)*dir.x - (x - x0)*dir.z]^2 + +//! [(x - x0)*dir.y - (y - y0)*dir.x]^2 +//! +//! First derivative of square vector product are: +//! Dx0 = 2*[(z - z0)*dir.x - (x - x0)*dir.z]*dir.z +//! -2*[(x - x0)*dir.y - (y - y0)*dir.x]*dir.y +//! Dy0 = -2*[(y - y0)*dir.z - (z - z0)*dir.y]*dir.z +//! +2*[(x - x0)*dir.y - (y - y0)*dir.x]*dir.x +//! Dz0 = 2*[(y - y0)*dir.z - (z - z0)*dir.y]*dir.y +//! -2*[(z - z0)*dir.x - (x - x0)*dir.z]*dir.x +//! +//! dF/dx0 : G1(...) = 2*Sum{[...]*Dx0} +//! dF/dy0 : G2(...) = 2*Sum{[...]*Dy0} +//! dF/dz0 : G3(...) = 2*Sum{[...]*Dz0} +//! dF/dR : G4(...) = -4*R*Sum[...] +//! [...] = [|(P(i) - Loc)^dir|^2 - R^2] +class ShapeAnalysis_FuncCylinderLSDist : public math_MultipleVarFunctionWithGradient +{ +public: + + DEFINE_STANDARD_ALLOC + + //! Constructor. + Standard_EXPORT ShapeAnalysis_FuncCylinderLSDist() {}; + + Standard_EXPORT ShapeAnalysis_FuncCylinderLSDist(const Handle(TColgp_HArray1OfXYZ)& thePoints, + const gp_Dir& theDir); + + void SetPoints(const Handle(TColgp_HArray1OfXYZ)& thePoints) + { + myPoints = thePoints; + } + + void SetDir(const gp_Dir& theDir) + { + myDir = theDir; + } + + //! Number of variables. + Standard_EXPORT Standard_Integer NbVariables() const Standard_OVERRIDE; + + //! Value. + Standard_EXPORT Standard_Boolean Value(const math_Vector& X,Standard_Real& F) Standard_OVERRIDE; + + //! Gradient. + Standard_EXPORT Standard_Boolean Gradient(const math_Vector& X,math_Vector& G) Standard_OVERRIDE; + + //! Value and gradient. + Standard_EXPORT Standard_Boolean Values(const math_Vector& X,Standard_Real& F,math_Vector& G) Standard_OVERRIDE; + +private: + + Handle(TColgp_HArray1OfXYZ) myPoints; + gp_Dir myDir; + +}; +#endif // _ShapeAnalysis_FuncCylinderLSDist_HeaderFile diff --git a/src/ShapeAnalysis/ShapeAnalysis_FuncSphereLSDist.cxx b/src/ShapeAnalysis/ShapeAnalysis_FuncSphereLSDist.cxx new file mode 100644 index 0000000000..7bf64af460 --- /dev/null +++ b/src/ShapeAnalysis/ShapeAnalysis_FuncSphereLSDist.cxx @@ -0,0 +1,115 @@ +// Created on: 2016-05-10 +// Created by: Alexander MALYSHEV +// Copyright (c) 1995-1999 Matra Datavision +// Copyright (c) 1999-2016 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + + +#include +#include +#include +#include + +//======================================================================= +//function : ShapeAnalysis_FuncSphereLSDist +//purpose : +//======================================================================= +ShapeAnalysis_FuncSphereLSDist::ShapeAnalysis_FuncSphereLSDist(const Handle(TColgp_HArray1OfXYZ)& thePoints): + myPoints(thePoints) +{ +} + +//======================================================================= +//function : NbVariables +//purpose : +//======================================================================= +Standard_Integer ShapeAnalysis_FuncSphereLSDist::NbVariables () const +{ + return 4; +} + +//======================================================================= +//function : Value +//purpose : +//======================================================================= +Standard_Boolean ShapeAnalysis_FuncSphereLSDist::Value(const math_Vector& X,Standard_Real& F) +{ + gp_XYZ aLoc(X(1), X(2), X(3)); + Standard_Real anR2 = X(4)*X(4); + + F = 0.; + Standard_Integer i; + for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i) + { + Standard_Real d = (myPoints->Value(i) - aLoc).SquareModulus() - anR2; + F += d * d; + } + + return Standard_True; +} + +//======================================================================= +//function : Gradient +//purpose : +//======================================================================= +Standard_Boolean ShapeAnalysis_FuncSphereLSDist::Gradient(const math_Vector& X,math_Vector& G) + +{ + gp_XYZ aLoc(X(1), X(2), X(3)); + Standard_Real anR = X(4), anR2 = anR * anR; + + G.Init(0.); + Standard_Integer i; + for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i) + { + gp_XYZ dLoc = myPoints->Value(i) - aLoc; + Standard_Real d = dLoc.SquareModulus() - anR2; + G(1) += d * dLoc.X(); + G(2) += d * dLoc.Y(); + G(3) += d * dLoc.Z(); + G(4) += d; + } + G *= -4; + G(4) *= anR; + + return Standard_True; +} + +//======================================================================= +//function : Values +//purpose : +//======================================================================= +Standard_Boolean ShapeAnalysis_FuncSphereLSDist::Values(const math_Vector& X,Standard_Real& F,math_Vector& G) +{ + gp_XYZ aLoc(X(1), X(2), X(3)); + Standard_Real anR = X(4), anR2 = anR * anR; + + G.Init(0.); + F = 0.; + Standard_Integer i; + for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i) + { + gp_XYZ dLoc = myPoints->Value(i) - aLoc; + Standard_Real d = dLoc.SquareModulus() - anR2; + G(1) += d * dLoc.X(); + G(2) += d * dLoc.Y(); + G(3) += d * dLoc.Z(); + G(4) += d; + F += d * d; + } + G *= -4; + G(4) *= anR; + + return true; +} + diff --git a/src/ShapeAnalysis/ShapeAnalysis_FuncSphereLSDist.hxx b/src/ShapeAnalysis/ShapeAnalysis_FuncSphereLSDist.hxx new file mode 100644 index 0000000000..fae0ffa8c6 --- /dev/null +++ b/src/ShapeAnalysis/ShapeAnalysis_FuncSphereLSDist.hxx @@ -0,0 +1,77 @@ + +// Copyright (c) 1991-1999 Matra Datavision +// Copyright (c) 1999-2022 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _ShapeAnalysis_FuncSphereLSDist_HeaderFile +#define _ShapeAnalysis_FuncSphereLSDist_HeaderFile + +#include +#include + +#include +#include +#include + +//! Function for search of sphere canonic parameters: coordinates of center and radius from set of moints +//! by least square method. +//! //! +//! The class inherits math_MultipleVarFunctionWithGradient and thus is intended +//! for use in math_BFGS algorithm. +//! +//! The criteria is: +//! F(x0, y0, z0, R) = Sum[(x(i) - x0)^2 + (y(i) - y0)^2 + (z(i) - z0)^2 - R^2]^2 => min, +//! x(i), y(i), z(i) - coordinates of sample points, x0, y0, z0, R - coordinates of center and radius of sphere, +//! which must be defined +//! +//! The first derivative are: +//! dF/dx0 : G1(x0, y0, z0, R) = -4*Sum{[...]*(x(i) - x0)} +//! dF/dy0 : G2(x0, y0, z0, R) = -4*Sum{[...]*(y(i) - y0)} +//! dF/dz0 : G3(x0, y0, z0, R) = -4*Sum{[...]*(z(i) - z0)} +//! dF/dR : G4(x0, y0, z0, R) = -4*R*Sum[...] +//! [...] = [(x(i) - x0)^2 + (y(i) - y0)^2 + (z(i) - z0)^2 - R^2] +//! +class ShapeAnalysis_FuncSphereLSDist : public math_MultipleVarFunctionWithGradient +{ +public: + + DEFINE_STANDARD_ALLOC + + //! Constructor. + Standard_EXPORT ShapeAnalysis_FuncSphereLSDist() {}; + + Standard_EXPORT ShapeAnalysis_FuncSphereLSDist(const Handle(TColgp_HArray1OfXYZ)& thePoints); + + void SetPoints(const Handle(TColgp_HArray1OfXYZ)& thePoints) + { + myPoints = thePoints; + } + + //! Number of variables. + Standard_EXPORT Standard_Integer NbVariables() const Standard_OVERRIDE; + + //! Value. + Standard_EXPORT Standard_Boolean Value(const math_Vector& X,Standard_Real& F) Standard_OVERRIDE; + + //! Gradient. + Standard_EXPORT Standard_Boolean Gradient(const math_Vector& X,math_Vector& G) Standard_OVERRIDE; + + //! Value and gradient. + Standard_EXPORT Standard_Boolean Values(const math_Vector& X,Standard_Real& F,math_Vector& G) Standard_OVERRIDE; + +private: + + Handle(TColgp_HArray1OfXYZ) myPoints; + +}; +#endif // _ShapeAnalysis_FuncSphereLSDist_HeaderFile diff --git a/src/gce/gce_MakeCirc.cxx b/src/gce/gce_MakeCirc.cxx index 3e74a49982..18ecf6b2b1 100644 --- a/src/gce/gce_MakeCirc.cxx +++ b/src/gce/gce_MakeCirc.cxx @@ -53,7 +53,7 @@ gce_MakeCirc::gce_MakeCirc(const gp_Pnt& P1 , //========================================================================= Standard_Real dist1, dist2, dist3, aResolution; // - aResolution=gp::Resolution(); + aResolution = gp::Resolution(); // dist1 = P1.Distance(P2); dist2 = P1.Distance(P3); @@ -76,7 +76,7 @@ gce_MakeCirc::gce_MakeCirc(const gp_Pnt& P1 , P2.Coord(x2,y2,z2); P3.Coord(x3,y3,z3); gp_Dir Dir1(x2-x1,y2-y1,z2-z1); - gp_Dir Dir2(x3-x2,y3-y2,z3-z2); + gp_Vec VDir2(x3-x2,y3-y2,z3-z2); // gp_Ax1 anAx1(P1, Dir1); gp_Lin aL12 (anAx1); @@ -85,14 +85,21 @@ gce_MakeCirc::gce_MakeCirc(const gp_Pnt& P1 , return; } // - gp_Dir Dir3 = Dir1.Crossed(Dir2); + gp_Vec VDir1(Dir1); + gp_Vec VDir3 = VDir1.Crossed(VDir2); + if(VDir3.SquareMagnitude() < aResolution) + { + TheError = gce_ColinearPoints; + return; + } // + gp_Dir Dir3(VDir3); gp_Dir dir = Dir1.Crossed(Dir3); gp_Lin L1(gp_Pnt((P1.XYZ()+P2.XYZ())/2.),dir); - dir = Dir2.Crossed(Dir3); + dir = VDir2.Crossed(Dir3); gp_Lin L2(gp_Pnt((P3.XYZ()+P2.XYZ())/2.),dir); - Standard_Real Tol = 0.000000001; + Standard_Real Tol = Precision::PConfusion(); Extrema_ExtElC distmin(L1,L2,Tol); if (!distmin.IsDone()) { @@ -139,7 +146,7 @@ gce_MakeCirc::gce_MakeCirc(const gp_Pnt& P1 , //Dir2 = gp_Dir(x2-x3,y2-y3,z2-z3); //modified by NIZNHY-PKV Thu Mar 3 11:31:13 2005t // - TheCirc = gp_Circ(gp_Ax2(pInt, Dir3, Dir1),(dist1+dist2+dist3)/3.); + TheCirc = gp_Circ(gp_Ax2(pInt, gp_Dir(VDir3), Dir1),(dist1+dist2+dist3)/3.); TheError = gce_Done; } } diff --git a/tests/cr/approx/A1 b/tests/cr/approx/A1 new file mode 100644 index 0000000000..f449dafb29 --- /dev/null +++ b/tests/cr/approx/A1 @@ -0,0 +1,12 @@ +polyline w 0 0 0 1 0 0 1 1 0 0 1 0 0 0 0 +getanasurf res w pln 1.e-3 +if {[isdraw res]} { + set log [dump res] + if { [regexp {Plane} $log ] != 1 } { + puts "Error: surface is not a plane" + } +} else { + puts "Error: required surface is not got" +} + + diff --git a/tests/cr/approx/A2 b/tests/cr/approx/A2 new file mode 100644 index 0000000000..5ac2f8f5c0 --- /dev/null +++ b/tests/cr/approx/A2 @@ -0,0 +1,15 @@ +cylinder cyl 0 0 0 1 1 1 1 +mkface f cyl 0 2 0 1 +explode f w +nurbsconvert w f_1 +cylinder cyl 0.1 0.1 0.1 1 1 1 .9 +getanasurf res w cyl 1.e-3 cyl +if {[isdraw res]} { + set log [dump res] + if { [regexp {CylindricalSurface} $log ] != 1 } { + puts "Error: surface is not a cylindrical surface" + } +} else { + puts "Error: required surface is not got" +} + diff --git a/tests/cr/approx/A3 b/tests/cr/approx/A3 new file mode 100644 index 0000000000..f59b3c6305 --- /dev/null +++ b/tests/cr/approx/A3 @@ -0,0 +1,15 @@ +plane pl 0 0 0 1 1 1 +pcone con pl 2 1 3 60 +explode con w +nurbsconvert w con_1 +cone con 0.1 .1 .1 1 1 1 -17 1.9 +getanasurf res w con 1.e-3 con +if {[isdraw res]} { + set log [dump res] + if { [regexp {ConicalSurface} $log ] != 1 } { + puts "Error: surface is not a conical surface" + } +} else { + puts "Error: required surface is not got" +} + diff --git a/tests/cr/approx/A4 b/tests/cr/approx/A4 new file mode 100644 index 0000000000..d1c3384f30 --- /dev/null +++ b/tests/cr/approx/A4 @@ -0,0 +1,16 @@ +sphere sph 1 +mkface f sph 0 2 0 pi/2 +explode f w +nurbsconvert w f_1 +sphere sph 0.1 0.1 0.1 1 +getanasurf res w sph 1.e-3 sph +if {[isdraw res]} { + set log [dump res] + if { [regexp {SphericalSurface} $log ] != 1 } { + puts "Error: surface is not a spherical surface" + } +} else { + puts "Error: required surface is not got" +} + + diff --git a/tests/cr/base/A1 b/tests/cr/base/A1 new file mode 100644 index 0000000000..4de2c8d7bb --- /dev/null +++ b/tests/cr/base/A1 @@ -0,0 +1,13 @@ +beziercurve bz 3 0 0 0 .5 .5e-3 0 1 0 0 +mkedge e bz +getanacurve c e lin 1.e-3 +if {[isdraw c]} { + set log [dump c] + if { [regexp {Line} $log ] != 1 } { + puts "Error: curve is not a line" + } +} else { + puts "Error: required curve is not got" +} + + diff --git a/tests/cr/base/A2 b/tests/cr/base/A2 new file mode 100644 index 0000000000..5f66a23ee4 --- /dev/null +++ b/tests/cr/base/A2 @@ -0,0 +1,12 @@ +beziercurve bz 3 0 0 0 .5 .5e-3 0 1 0 0 +mkedge e bz +getanacurve c e cir 1.e-3 +if {[isdraw c]} { + set log [dump c] + if { [regexp {Circle} $log ] != 1 } { + puts "Error: curve is not a circle" + } +} else { + puts "Error: required curve is not got" +} + diff --git a/tests/cr/base/A3 b/tests/cr/base/A3 new file mode 100644 index 0000000000..97d6e1ec5d --- /dev/null +++ b/tests/cr/base/A3 @@ -0,0 +1,13 @@ +ellipse bz 0 0 0 0 0 1 1 .5 +convert bz bz +mkedge e bz +getanacurve c e ell 1.e-7 +if {[isdraw c]} { + set log [dump c] + if { [regexp {Ellipse} $log ] != 1 } { + puts "Error: curve is not an ellipse" + } +} else { + puts "Error: required curve is not got" +} + diff --git a/tests/cr/base/A4 b/tests/cr/base/A4 new file mode 100644 index 0000000000..76c26f34b4 --- /dev/null +++ b/tests/cr/base/A4 @@ -0,0 +1,15 @@ +beziercurve bz 3 0 0 0 .5 .5e-6 0 1 0 0 +mkedge e1 bz 0 .3 +mkedge e2 bz .3 0.7 +mkedge e3 bz .7 1. +wire w e1 e2 e3 +getanacurve c w lin 1.e-3 +if {[isdraw c]} { + set log [dump c] + if { [regexp {Line} $log ] != 1 } { + puts "Error: curve is not a line" + } +} else { + puts "Error: required curve is not got" +} + diff --git a/tests/cr/base/A5 b/tests/cr/base/A5 new file mode 100644 index 0000000000..7323a6221a --- /dev/null +++ b/tests/cr/base/A5 @@ -0,0 +1,15 @@ +circle bz 0 0 0 0 0 1 1 +convert bz bz +mkedge e1 bz 0 1 +mkedge e2 bz 1 2.5 +mkedge e3 bz 2.5 6 +wire w e1 e2 e3 +getanacurve c w cir 1.e-7 +if {[isdraw c]} { + set log [dump c] + if { [regexp {Circle} $log ] != 1 } { + puts "Error: curve is not a circle" + } +} else { + puts "Error: required curve is not got" +} diff --git a/tests/cr/base/A6 b/tests/cr/base/A6 new file mode 100644 index 0000000000..dabfd08c3e --- /dev/null +++ b/tests/cr/base/A6 @@ -0,0 +1,15 @@ +ellipse bz 0 0 0 0 0 1 1 .5 +convert bz bz +mkedge e1 bz 0 1 +mkedge e2 bz 1 2.5 +mkedge e3 bz 2.5 6 +wire w e1 e2 e3 +getanacurve c w ell 1.e-7 +if {[isdraw c]} { + set log [dump c] + if { [regexp {Ellipse} $log ] != 1 } { + puts "Error: curve is not an ellipse" + } +} else { + puts "Error: required curve is not got" +} diff --git a/tests/cr/base/B1 b/tests/cr/base/B1 new file mode 100644 index 0000000000..87498c8809 --- /dev/null +++ b/tests/cr/base/B1 @@ -0,0 +1,15 @@ +plane surf +trim surf surf -1 1 -1 1 +convert surf surf +mkface f surf +getanasurf asurf f pln 1.e-7 +if {[isdraw asurf]} { + set log [dump asurf] + if { [regexp {Plane} $log ] != 1 } { + puts "Error: surface is not a plane" + } +} else { + puts "Error: required surface is not got" +} + + diff --git a/tests/cr/base/B2 b/tests/cr/base/B2 new file mode 100644 index 0000000000..312e25ab4b --- /dev/null +++ b/tests/cr/base/B2 @@ -0,0 +1,14 @@ +cylinder surf 1 +trimv surf surf -1 1 +convert surf surf +mkface f surf +getanasurf asurf f cyl 1.e-7 +if {[isdraw asurf]} { + set log [dump asurf] + if { [regexp {CylindricalSurface} $log ] != 1 } { + puts "Error: surface is not a cylindrical surface" + } +} else { + puts "Error: required surface is not got" +} + diff --git a/tests/cr/base/B3 b/tests/cr/base/B3 new file mode 100644 index 0000000000..ce1e30b4b1 --- /dev/null +++ b/tests/cr/base/B3 @@ -0,0 +1,13 @@ +cone surf 30 0 +trimv surf surf -1 0 +convert surf surf +mkface f surf +getanasurf asurf f con 1.e-7 +if {[isdraw asurf]} { + set log [dump asurf] + if { [regexp {ConicalSurface} $log ] != 1 } { + puts "Error: surface is not a conical surface" + } +} else { + puts "Error: required surface is not got" +} diff --git a/tests/cr/base/B4 b/tests/cr/base/B4 new file mode 100644 index 0000000000..7259721402 --- /dev/null +++ b/tests/cr/base/B4 @@ -0,0 +1,12 @@ +sphere surf 1 +convert surf surf +mkface f surf +getanasurf asurf f sph 1.e-7 +if {[isdraw asurf]} { + set log [dump asurf] + if { [regexp {SphericalSurface} $log ] != 1 } { + puts "Error: surface is not a spherical surface" + } +} else { + puts "Error: required surface is not got" +} diff --git a/tests/cr/base/B5 b/tests/cr/base/B5 new file mode 100644 index 0000000000..58c4d5b223 --- /dev/null +++ b/tests/cr/base/B5 @@ -0,0 +1,25 @@ +plane surf +trim surf1 surf -1 0 -1 0 +convert surf1 surf1 +mkface f1 surf1 +trim surf2 surf -1 0 0 1 +convert surf2 surf2 +mkface f2 surf2 +trim surf3 surf 0 1 0 1 +convert surf3 surf3 +mkface f3 surf3 +trim surf4 surf 0 1 -1 0 +convert surf4 surf4 +mkface f4 surf4 +sewing sh f1 f2 f3 f4 +getanasurf asurf sh pln 1.e-7 +if {[isdraw asurf]} { + set log [dump asurf] + if { [regexp {Plane} $log ] != 1 } { + puts "Error: surface is not a plane" + } +} else { + puts "Error: required surface is not got" +} + + diff --git a/tests/cr/base/B6 b/tests/cr/base/B6 new file mode 100644 index 0000000000..db670f3a5d --- /dev/null +++ b/tests/cr/base/B6 @@ -0,0 +1,25 @@ +cylinder surf 1 +trim surf1 surf 0 3 -1 0 +convert surf1 surf1 +mkface f1 surf1 +trim surf2 surf 0 3 0 1 +convert surf2 surf2 +mkface f2 surf2 +trim surf3 surf 3 6 0 1 +convert surf3 surf3 +mkface f3 surf3 +trim surf4 surf 3 6 -1 0 +convert surf4 surf4 +mkface f4 surf4 +sewing sh f1 f2 f3 f4 +getanasurf asurf sh cyl 1.e-7 +if {[isdraw asurf]} { + set log [dump asurf] + if { [regexp {CylindricalSurface} $log ] != 1 } { + puts "Error: surface is not a cylindrical surface" + } +} else { + puts "Error: required surface is not got" +} + + diff --git a/tests/cr/base/B7 b/tests/cr/base/B7 new file mode 100644 index 0000000000..2c2d5f2fc4 --- /dev/null +++ b/tests/cr/base/B7 @@ -0,0 +1,25 @@ +cone surf 30 0 +trim surf1 surf 0 3 0 1 +convert surf1 surf1 +mkface f1 surf1 +trim surf2 surf 0 3 1 2 +convert surf2 surf2 +mkface f2 surf2 +trim surf3 surf 3 6 1 2 +convert surf3 surf3 +mkface f3 surf3 +trim surf4 surf 3 6 0 1 +convert surf4 surf4 +mkface f4 surf4 +sewing sh f1 f2 f3 f4 +getanasurf asurf sh con 1.e-7 +if {[isdraw asurf]} { + set log [dump asurf] + if { [regexp {ConicalSurface} $log ] != 1 } { + puts "Error: surface is not a conical surface" + } +} else { + puts "Error: required surface is not got" +} + + diff --git a/tests/cr/base/B8 b/tests/cr/base/B8 new file mode 100644 index 0000000000..64e69c9959 --- /dev/null +++ b/tests/cr/base/B8 @@ -0,0 +1,25 @@ +sphere surf 1 +trim surf1 surf 0 3 -1.5 0 +convert surf1 surf1 +mkface f1 surf1 +trim surf2 surf 0 3 0 1.5 +convert surf2 surf2 +mkface f2 surf2 +trim surf3 surf 3 6 0 1.5 +convert surf3 surf3 +mkface f3 surf3 +trim surf4 surf 3 6 -1.5 0 +convert surf4 surf4 +mkface f4 surf4 +sewing sh f1 f2 f3 f4 +getanasurf asurf sh sph 1.e-7 +if {[isdraw asurf]} { + set log [dump asurf] + if { [regexp {SphericalSurface} $log ] != 1 } { + puts "Error: surface is not a spherical surface" + } +} else { + puts "Error: required surface is not got" +} + + diff --git a/tests/cr/base/B9 b/tests/cr/base/B9 new file mode 100644 index 0000000000..195c30074a --- /dev/null +++ b/tests/cr/base/B9 @@ -0,0 +1,40 @@ +cylinder surf 1 +trim surf1 surf 0 1 -1 1 +convert surf1 surf1 +mkface f1 surf1 +trim surf2 surf 1 2 -1 1 +convert surf2 surf2 +mkface f2 surf2 +trim surf3 surf 2 3 -1 1 +convert surf3 surf3 +mkface f3 surf3 +trim surf4 surf 3 4 -1 1 +convert surf4 surf4 +mkface f4 surf4 +sewing sh f1 f2 f3 f4 +cylinder cyl 0 0 -2.5 1 0 0 3 +trimv cyl1 cyl -2 2 +mkface fcyl cyl1 +nurbsconvert fcyl fcyl +bsection ss fcyl sh +explode ss e +shape w w +add ss_1 w +add ss_2 w +add ss_3 w +add ss_4 w +add ss_5 w +getanasurf asurf w cyl 1.e-7 cyl +if {[isdraw asurf]} { + set log [dump asurf] + if { [regexp {CylindricalSurface} $log ] != 1 } { + puts "Error: surface is not a cylindrical surface" + } + if { [regexp {Radius :3} $log ] != 1 } { + puts "Error: surface is not a sample surface" + } +} else { + puts "Error: required surface is not got" +} + + diff --git a/tests/cr/begin b/tests/cr/begin new file mode 100644 index 0000000000..f86d1068c6 --- /dev/null +++ b/tests/cr/begin @@ -0,0 +1,14 @@ +if { [array get Draw_Groups "TOPOLOGY Check commands"] == "" } { + pload TOPTEST +} + +# To prevent loops limit to 5 minutes +cpulimit 300 + +if { [info exists imagedir] == 0 } { + set imagedir . +} + +if { [info exists test_image ] == 0 } { + set test_image photo +} diff --git a/tests/cr/end b/tests/cr/end new file mode 100644 index 0000000000..601b1e1213 --- /dev/null +++ b/tests/cr/end @@ -0,0 +1,3 @@ + +# to end a test script +puts "TEST COMPLETED" diff --git a/tests/cr/grids.list b/tests/cr/grids.list new file mode 100644 index 0000000000..8eb5c913e5 --- /dev/null +++ b/tests/cr/grids.list @@ -0,0 +1,2 @@ +001 base +002 approx