1
0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-08-19 13:40:49 +03:00

0032995: Modeling algorithms - Interface for checking canonical geometry

GeomConvert_CurveToAnaCurve, GeomConvert_SurfToAnaSurf - geometrical algorithms
for converting geometrical curve and surfaces in canonical geometry with given tolerance.

ShapeAnalysis_CanonicalRecognition - interface for checking canonical geometry.

# Conflicts:
#	src/GeomConvert/FILES
#	src/GeomConvert/GeomConvert_CurveToAnaCurve.cxx
#	src/GeomConvert/GeomConvert_CurveToAnaCurve.hxx
#	src/GeomConvert/GeomConvert_SurfToAnaSurf.cxx
#	src/GeomConvert/GeomConvert_SurfToAnaSurf.hxx
#	src/SWDRAW/SWDRAW_ShapeAnalysis.cxx
#	src/ShapeAnalysis/FILES
#	src/ShapeAnalysis/ShapeAnalysis_CanonicalRecognition.cxx
#	src/ShapeAnalysis/ShapeAnalysis_CanonicalRecognition.hxx
#	tests/cr/base/A1
#	tests/cr/base/A2
#	tests/cr/base/A3
#	tests/cr/base/A4
#	tests/cr/base/A5
#	tests/cr/base/A6
#	tests/cr/base/B1
#	tests/cr/base/B2
#	tests/cr/base/B3
#	tests/cr/base/B4
#	tests/cr/base/B5
#	tests/cr/base/B6
#	tests/cr/base/B7
#	tests/cr/base/B8
#	tests/cr/base/B9
#	tests/cr/grids.list
This commit is contained in:
vsv
2022-07-07 16:39:51 +03:00
parent 2cdc84e120
commit dc21bf35fe
37 changed files with 1846 additions and 800 deletions

View File

@@ -896,6 +896,25 @@ Analysis of surfaces is allowed for following shapes:
* 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.
@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. This means that not only it corrects and upgrades but also changes the definition of a shape with regard to its geometry, size and other aspects. Convenient API allows you to create your own tools to perform specific upgrading. Additional tools for particular cases provide an ability to divide shapes and surfaces according to certain criteria.

View File

@@ -37,34 +37,10 @@
#include <TColStd_Array1OfReal.hxx>
#include <TColStd_Array2OfReal.hxx>
#include <GeomAbs_CurveType.hxx>
#include <math_Vector.hxx>
#include <math_Matrix.hxx>
#include <math_Gauss.hxx>
//=======================================================================
//function : AdjustByPeriod
//purpose :
//=======================================================================
static Standard_Real AdjustByPeriod(const Standard_Real Val,
const Standard_Real ToVal,
const Standard_Real Period)
{
Standard_Real diff = Val - ToVal;
Standard_Real D = Abs(diff);
Standard_Real P = Abs(Period);
if (D <= 0.5 * P) return 0.;
if (P < 1e-100) return diff;
return (diff >0 ? -P : P) * floor(D / P + 0.5);
}
//=======================================================================
//function : AdjustToPeriod
//purpose :
//=======================================================================
static Standard_Real AdjustToPeriod(const Standard_Real Val,
const Standard_Real ValMin,
const Standard_Real ValMax)
{
return AdjustByPeriod(Val, 0.5 * (ValMin + ValMax), ValMax - ValMin);
}
GeomConvert_CurveToAnaCurve::GeomConvert_CurveToAnaCurve():
myGap(Precision::Infinite()),
@@ -141,8 +117,7 @@ Standard_Boolean GeomConvert_CurveToAnaCurve::IsLinear(const TColgp_Array1OfPnt&
}
}
Standard_Real dPreci = Precision::Confusion()*Precision::Confusion();
if(dMax < dPreci)
if (dMax < Precision::SquareConfusion())
return Standard_False;
Standard_Real tol2 = tolerance*tolerance;
@@ -197,8 +172,7 @@ Handle(Geom_Line) GeomConvert_CurveToAnaCurve::ComputeLine (const Handle(Geom_Cu
gp_Pnt P1 = curve->Value (c1);
gp_Pnt P2 = curve->Value (c2);
Standard_Real dPreci = Precision::Confusion()*Precision::Confusion();
if(P1.SquareDistance(P2) < dPreci)
if(P1.SquareDistance(P2) < Precision::SquareConfusion())
return line;
cf = c1; cl = c2;
@@ -244,11 +218,9 @@ Handle(Geom_Line) GeomConvert_CurveToAnaCurve::ComputeLine (const Handle(Geom_Cu
//=======================================================================
Standard_Boolean GeomConvert_CurveToAnaCurve::GetCircle (gp_Circ& crc,
const gp_Pnt& P0,const gp_Pnt& P1, const gp_Pnt& P2,
const Standard_Real d0, const Standard_Real d1, const Standard_Real epsang)
const gp_Pnt& P0,const gp_Pnt& P1, const gp_Pnt& P2)
{
// Control if points are not aligned (should be done by MakeCirc ?)
// d0 and d1 are distances (p0 p1) and (p0 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;
@@ -256,15 +228,12 @@ Standard_Boolean GeomConvert_CurveToAnaCurve::GetCircle (gp_Circ& crc,
return Standard_False;
if (Abs(P2.X()) > aMaxCoord || Abs(P2.Y()) > aMaxCoord || Abs(P2.Z()) > aMaxCoord)
return Standard_False;
gp_Vec p0p1 (P0,P1);
gp_Vec p0p2 (P0,P2);
Standard_Real ang = p0p1.CrossSquareMagnitude (p0p2);
if (ang < d0*d1*epsang) 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();
@@ -299,22 +268,9 @@ Handle(Geom_Curve) GeomConvert_CurveToAnaCurve::ComputeCircle (const Handle(Geom
P0 = c3d->Value(c1);
P1 = c3d->Value(ca);
P2 = c3d->Value(cb);
// For Control if points are not aligned (should be done by MakeCirc ?),
// d0 and d1 are distances (p0 p1) and (p0 p2)
// return result here means : void result (no circle)
Standard_Real eps = Precision::Angular(); // angular resolution
Standard_Real d0 = P0.Distance(P1); Standard_Real d1 = P0.Distance(P2);
if (d0 < Precision::Confusion() || d1 < Precision::Confusion()) return circ;
gp_Circ crc;
if (!GetCircle (crc,P0,P1,P2,d0,d1,eps)) return circ;
cf = 0;
// Standard_Real cm = ElCLib::Parameter (crc,c3d->Value ((c1+c2)/2) );
// cl = ElCLib::Parameter (crc,c3d->Value (c2));
// if (cl > 2 * M_PI) cl = 2 * M_PI;
// if (cl < cm && cm > 0.99 * M_PI) cl = 2 * M_PI;
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;
@@ -334,17 +290,19 @@ Handle(Geom_Curve) GeomConvert_CurveToAnaCurve::ComputeCircle (const Handle(Geom
Standard_Real PI2 = 2 * M_PI;
cf = ElCLib::Parameter (crc,c3d->Value (c1));
cf+= AdjustToPeriod(cf,0.,PI2);
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+= AdjustToPeriod(cm,cf,cf+PI2);
cm = ElCLib::InPeriod(cm, cf, cf + PI2);
cl = ElCLib::Parameter (crc,c3d->Value (c2));
cl+= AdjustToPeriod(cl,cm,cm+PI2);
cl = ElCLib::InPeriod(cl, cm, cm + PI2);
circ = new Geom_Circle (crc);
return circ;
}
@@ -404,117 +362,6 @@ static Standard_Boolean IsArrayPntPlanar(const Handle(TColgp_HArray1OfPnt)& HAP,
return Standard_True;
}
//=======================================================================
//function : Determ3
//purpose :
//=======================================================================
static Standard_Real Determ3(const TColStd_Array2OfReal& A)
{
Standard_Real det = A.Value(1,1)*A.Value(2,2)*A.Value(3,3) +
A.Value(1,2)*A.Value(2,3)*A.Value(3,1) +
A.Value(1,3)*A.Value(2,1)*A.Value(3,2) -
A.Value(1,3)*A.Value(2,2)*A.Value(3,1) -
A.Value(1,1)*A.Value(2,3)*A.Value(3,2) -
A.Value(1,2)*A.Value(2,1)*A.Value(3,3);
return det;
}
//=======================================================================
//function : Determ4
//purpose :
//=======================================================================
static Standard_Real Determ4(const TColStd_Array2OfReal& A)
{
Standard_Real det=0;
Standard_Integer i,j,k;
TColStd_Array2OfReal C(1,3,1,3);
for(j=1; j<=4; j++) {
for(i=2; i<=4; i++) {
for(k=1; k<=4; k++) {
if(k<j)
C.SetValue(i-1,k,A.Value(i,k));
if(k>j)
C.SetValue(i-1,k-1,A.Value(i,k));
}
}
if(j==1 || j==3)
det = det + A.Value(1,j)*Determ3(C);
else
det = det - A.Value(1,j)*Determ3(C);
}
return det;
}
//=======================================================================
//function : Determ5
//purpose :
//=======================================================================
static Standard_Real Determ5(const TColStd_Array2OfReal& A)
{
Standard_Real det=0;
Standard_Integer i,j,k;
TColStd_Array2OfReal C(1,4,1,4);
for(j=1; j<=5; j++) {
for(i=2; i<=5; i++) {
for(k=1; k<=5; k++) {
if(k<j)
C.SetValue(i-1,k,A.Value(i,k));
if(k>j)
C.SetValue(i-1,k-1,A.Value(i,k));
}
}
if(j==1 || j==3 || j==5)
det = det + A.Value(1,j)*Determ4(C);
else
det = det - A.Value(1,j)*Determ4(C);
}
return det;
}
//=======================================================================
//function : SolveLinSys5
//purpose :
//=======================================================================
static Standard_Boolean SolveLinSys5(const TColStd_Array2OfReal& A,
const TColStd_Array1OfReal& B,
TColStd_Array1OfReal& X )
{
Standard_Integer i,j;
Standard_Real D0,DN;
TColStd_Array2OfReal A1(1,5,1,5);
TColStd_Array1OfReal Tmp(1,5);
for(j=1; j<=5; j++) {
for(i=1; i<=5; i++) {
A1.SetValue(i,j,A.Value(i,j));
}
}
D0 = Determ5(A1);
if(D0==0)
return Standard_False;
for(j=1; j<=5; j++) {
for(i=1; i<=5; i++) {
Tmp.SetValue(i,A1.Value(i,j));
A1.SetValue(i,j,B.Value(i));
}
DN = Determ5(A1);
X.SetValue(j,DN/D0);
for(i=1; i<=5; i++) {
A1.SetValue(i,j,Tmp.Value(i));
}
}
return Standard_True;
}
//=======================================================================
//function : ConicdDefinition
//purpose :
@@ -675,30 +522,33 @@ Handle(Geom_Curve) GeomConvert_CurveToAnaCurve::ComputeEllipse(const Handle(Geom
Tr.SetTransformation(AX);
gp_Trsf Tr2 = Tr.Inverted();
TColStd_Array2OfReal Dt(1,5,1,5);
TColStd_Array1OfReal F(1,5), Sl(1,5);
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.SetValue(i, 1, XN*XN);
Dt.SetValue(i,2,XN*YN);
Dt.SetValue(i,3,YN*YN);
Dt.SetValue(i,4,XN);
Dt.SetValue(i,5,YN);
F.SetValue(i,-1);
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.;
}
//
if(!SolveLinSys5(Dt,F,Sl))
math_Gauss aSolver(Dt);
if (!aSolver.IsDone())
return res;
aSolver.Solve(F, Sl);
AF=Sl.Value(1);
BF=Sl.Value(2);
CF=Sl.Value(3);
DF=Sl.Value(4);
EF=Sl.Value(5);
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;
@@ -752,21 +602,23 @@ if (Q2 > 0 && Q1*Q3 < 0) {
Standard_Real PI2 = 2 * M_PI;
cf = ElCLib::Parameter(anEllipse, c3d->Value(c1));
cf += AdjustToPeriod(cf, 0., PI2);
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 += AdjustToPeriod(cm, cf, cf + PI2);
cm = ElCLib::InPeriod(cm, cf, cf + PI2);
cl = ElCLib::Parameter(anEllipse, c3d->Value(c2));
cl += AdjustToPeriod(cl, cm, cm + PI2);
cl = ElCLib::InPeriod(cl, cm, cm + PI2);
res = gell;
}
}
/*
if (Q2 < 0 && Q1 != 0) {
// hyberbola
}
@@ -774,7 +626,7 @@ if (Q2 < 0 && Q1 != 0) {
if (Q2 == 0 && Q1 != 0) {
// parabola
}
*/
return res;
}

View File

@@ -79,7 +79,7 @@ public:
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, const Standard_Real d0, const Standard_Real d1, const Standard_Real eps);
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

View File

@@ -30,9 +30,10 @@
#include <Geom_Plane.hxx>
#include <Geom_SphericalSurface.hxx>
#include <Geom_Surface.hxx>
#include <Geom_SurfaceOfRevolution.hxx>
#include <Geom_ToroidalSurface.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <GeomAdaptor_HSurface.hxx>
#include <GeomAdaptor_Surface.hxx>
#include <GeomAdaptor_Surface.hxx>
#include <GeomLib_IsPlanarSurface.hxx>
#include <gp_Ax3.hxx>
@@ -45,7 +46,56 @@
#include <GeomConvert_SurfToAnaSurf.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <Extrema_ExtElC.hxx>
//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 <isTryUMajor> = Standard_True try to use V isoline radius as minor radaius.
@@ -71,12 +121,12 @@ static Handle(Geom_Surface) TryTorusSphere(const Handle(Geom_Surface)& theSurf,
if (isTryUMajor)
{
IsoCrv1 = theSurf->VIso(Param1 + ((Param2 - Param1) / 3.));
IsoCrv2 = theSurf->VIso(Param1 + ((Param2 - Param1)*2. / 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));
IsoCrv2 = theSurf->UIso(Param1 + ((Param2 - Param1) * 2. / 3));
}
Handle(Geom_Curve) Crv1 = GeomConvert_CurveToAnaCurve::ComputeCurve(IsoCrv1, toler, aParam1ToCrv, aParam2ToCrv, cf, cl,
@@ -103,7 +153,7 @@ static Handle(Geom_Surface) TryTorusSphere(const Handle(Geom_Surface)& theSurf,
aPnt2 = aCircle1->Circ().Location();
aPnt3 = aCircle2->Circ().Location();
Standard_Real eps = 1.e-09; // angular resolution
//Standard_Real eps = 1.e-09; // angular resolution
Standard_Real d0 = aPnt1.Distance(aPnt2); Standard_Real d1 = aPnt1.Distance(aPnt3);
gp_Circ circ;
@@ -120,7 +170,7 @@ static Handle(Geom_Surface) TryTorusSphere(const Handle(Geom_Surface)& theSurf,
return newSurface;
}
if (!GeomConvert_CurveToAnaCurve::GetCircle(circ, aPnt1, aPnt2, aPnt3, d0, d1, eps))
if (!GeomConvert_CurveToAnaCurve::GetCircle(circ, aPnt1, aPnt2, aPnt3) /*, d0, d1, eps)*/)
return newSurface;
Standard_Real aMajorR = circ.Radius();
@@ -131,7 +181,7 @@ static Handle(Geom_Surface) TryTorusSphere(const Handle(Geom_Surface)& theSurf,
return newSurface;
}
static Standard_Real ComputeGap(const Handle(Geom_Surface)& theSurf,
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())
{
@@ -142,26 +192,26 @@ static Standard_Real ComputeGap(const Handle(Geom_Surface)& theSurf,
gp_Cone aCon;
gp_Sphere aSphere;
gp_Torus aTor;
switch (aSType)
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;
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.;
@@ -177,9 +227,9 @@ static Standard_Real ComputeGap(const Handle(Geom_Surface)& theSurf,
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;
Standard_Real V = theV1 + DV * (j - 1) + DV2;
for (i = 1; i <= NP; i++) {
Standard_Real U = theU1 + DU*(i - 1) + DU2;
Standard_Real U = theU1 + DU * (i - 1) + DU2;
theSurf->D0(U, V, P3d);
switch (aSType) {
@@ -236,9 +286,9 @@ static Standard_Real ComputeGap(const Handle(Geom_Surface)& theSurf,
//purpose :
//=======================================================================
GeomConvert_SurfToAnaSurf::GeomConvert_SurfToAnaSurf() :
myGap (-1.),
myConvType (GeomConvert_Simplest),
myTarget (GeomAbs_Plane)
myGap(-1.),
myConvType(GeomConvert_Simplest),
myTarget(GeomAbs_Plane)
{
}
@@ -247,12 +297,12 @@ GeomConvert_SurfToAnaSurf::GeomConvert_SurfToAnaSurf() :
//purpose :
//=======================================================================
GeomConvert_SurfToAnaSurf::GeomConvert_SurfToAnaSurf (const Handle(Geom_Surface)& S) :
GeomConvert_SurfToAnaSurf::GeomConvert_SurfToAnaSurf(const Handle(Geom_Surface)& S) :
myGap(-1.),
myConvType(GeomConvert_Simplest),
myTarget(GeomAbs_Plane)
{
Init ( S );
Init(S);
}
//=======================================================================
@@ -260,7 +310,7 @@ GeomConvert_SurfToAnaSurf::GeomConvert_SurfToAnaSurf (const Handle(Geom_Surface)
//purpose :
//=======================================================================
void GeomConvert_SurfToAnaSurf::Init (const Handle(Geom_Surface)& S)
void GeomConvert_SurfToAnaSurf::Init(const Handle(Geom_Surface)& S)
{
mySurf = S;
}
@@ -270,19 +320,24 @@ void GeomConvert_SurfToAnaSurf::Init (const Handle(Geom_Surface)& S)
//purpose :
//=======================================================================
Handle(Geom_Surface) GeomConvert_SurfToAnaSurf::ConvertToAnalytical (const Standard_Real InitialToler)
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) ) {
if (Precision::IsInfinite(U1) && Precision::IsInfinite(U2)) {
U1 = -1.;
U2 = 1.;
}
if(Precision::IsInfinite(V1) && Precision::IsInfinite(V2) ) {
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);
return ConvertToAnalytical(InitialToler, U1, U2, V1, V2);
}
//=======================================================================
@@ -298,33 +353,33 @@ Handle(Geom_Surface) GeomConvert_SurfToAnaSurf::ConvertToAnalytical(const Standa
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;
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;
@@ -341,7 +396,7 @@ Handle(Geom_Surface) GeomConvert_SurfToAnaSurf::ConvertToAnalytical(const Standa
{
return newSurf[isurf];
}
else
else
{
if (Umin - U1 > aTolBnd)
{
@@ -363,7 +418,7 @@ Handle(Geom_Surface) GeomConvert_SurfToAnaSurf::ConvertToAnalytical(const Standa
V2 = Vmax;
aDoSegment = Standard_True;
}
}
Standard_Boolean IsBz = aSType == GeomAbs_BezierSurface;
@@ -399,8 +454,8 @@ Handle(Geom_Surface) GeomConvert_SurfToAnaSurf::ConvertToAnalytical(const Standa
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))
if (myConvType == GeomConvert_Simplest ||
(myConvType == GeomConvert_Target && myTarget == GeomAbs_Plane))
{
myGap = dd[isurf];
return newSurf[isurf];
@@ -424,10 +479,10 @@ Handle(Geom_Surface) GeomConvert_SurfToAnaSurf::ConvertToAnalytical(const Standa
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);
Standard_Real VMid = 0.5 * (V1 + V2);
Standard_Real UMid = 0.5 * (U1 + U2);
Handle(Geom_Surface) TrSurf = aTempS;
if(!aDoSegment)
if (!aDoSegment)
TrSurf = new Geom_RectangularTrimmedSurface(aTempS, U1, U2, V1, V2);
Handle(Geom_Curve) UIso = TrSurf->UIso(UMid);
@@ -566,7 +621,7 @@ Handle(Geom_Surface) GeomConvert_SurfToAnaSurf::ConvertToAnalytical(const Standa
gp_Ax3 Axes;
Standard_Real semiangle =
gp_Vec(isoline).Angle(gp_Vec(P3, P1));
if (semiangle>M_PI / 2) semiangle = M_PI - semiangle;
if (semiangle > M_PI / 2) semiangle = M_PI - semiangle;
if (R1 > R3) {
radius = R3;
Axes = gp_Ax3(P3, gp_Dir(gp_Vec(P3, P1)));
@@ -626,7 +681,7 @@ Handle(Geom_Surface) GeomConvert_SurfToAnaSurf::ConvertToAnalytical(const Standa
myGap = dd[imin];
return newSurf[imin];
}
return newSurf[0];
}
@@ -636,19 +691,19 @@ Handle(Geom_Surface) GeomConvert_SurfToAnaSurf::ConvertToAnalytical(const Standa
//=======================================================================
Standard_Boolean GeomConvert_SurfToAnaSurf::IsSame(const Handle(Geom_Surface)& S1,
const Handle(Geom_Surface)& S2,
const Standard_Real tol)
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)))
if (!S1->IsKind(STANDARD_TYPE(Geom_ElementarySurface)) ||
!S2->IsKind(STANDARD_TYPE(Geom_ElementarySurface)))
return Standard_False;
Handle(GeomAdaptor_HSurface) anAdaptor1 = new GeomAdaptor_HSurface(S1);
Handle(GeomAdaptor_HSurface) anAdaptor2 = new GeomAdaptor_HSurface(S2);
GeomAbs_SurfaceType aST1 = anAdaptor1->GetType();
GeomAbs_SurfaceType aST2 = anAdaptor2->GetType();
GeomAdaptor_Surface anAdaptor1(S1);
GeomAdaptor_Surface anAdaptor2(S2);
GeomAbs_SurfaceType aST1 = anAdaptor1.GetType();
GeomAbs_SurfaceType aST2 = anAdaptor2.GetType();
if (aST1 != aST2)
{
@@ -658,28 +713,28 @@ Standard_Boolean GeomConvert_SurfToAnaSurf::IsSame(const Handle(Geom_Surface)& S
IntAna_QuadQuadGeo interii;
if (aST1 == GeomAbs_Plane)
{
interii.Perform(anAdaptor1->Plane(), anAdaptor2->Plane(), tol, tol);
interii.Perform(anAdaptor1.Plane(), anAdaptor2.Plane(), tol, tol);
}
else if (aST1 == GeomAbs_Cylinder)
{
interii.Perform(anAdaptor1->Cylinder(), anAdaptor2->Cylinder(), tol);
interii.Perform(anAdaptor1.Cylinder(), anAdaptor2.Cylinder(), tol);
}
else if (aST1 == GeomAbs_Cone)
{
interii.Perform(anAdaptor1->Cone(), anAdaptor2->Cone(), tol);
interii.Perform(anAdaptor1.Cone(), anAdaptor2.Cone(), tol);
}
else if (aST1 == GeomAbs_Sphere)
{
interii.Perform(anAdaptor1->Sphere(), anAdaptor2->Sphere(), tol);
interii.Perform(anAdaptor1.Sphere(), anAdaptor2.Sphere(), tol);
}
else if (aST1 == GeomAbs_Torus)
{
interii.Perform(anAdaptor1->Torus(), anAdaptor2->Torus(), tol);
interii.Perform(anAdaptor1.Torus(), anAdaptor2.Torus(), tol);
}
if (!interii.IsDone())
return Standard_False;
IntAna_ResultType aTypeRes = interii.TypeInter();
return aTypeRes == IntAna_Same;
@@ -694,16 +749,16 @@ Standard_Boolean GeomConvert_SurfToAnaSurf::IsSame(const Handle(Geom_Surface)& S
Standard_Boolean GeomConvert_SurfToAnaSurf::IsCanonical(const Handle(Geom_Surface)& S)
{
if(S.IsNull())
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)))
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;
}

View File

@@ -62,14 +62,8 @@ public:
//! Tries to convert the Surface to an Analytic form
//! Returns the result
//! Works only if the Surface is BSpline or Bezier.
//! Else, or in case of failure, returns a Null Handle
//! In case of failure, returns a Null Handle
//!
//! If <substitute> is True, the new surface replaces the actual
//! one in <me>
//!
//! It works by analysing the case which can apply, creating the
//! corresponding analytic surface, then checking coincidence
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,

File diff suppressed because it is too large Load Diff

View File

@@ -994,7 +994,15 @@ static Standard_Integer getanasurf(Draw_Interpretor& di,
Standard_Integer n, const char** a)
{
if (n < 3) return 1;
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();
@@ -1007,8 +1015,17 @@ static Standard_Integer getanasurf(Draw_Interpretor& di,
GeomAbs_SurfaceType aTargets[] = { GeomAbs_Plane, GeomAbs_Cylinder, GeomAbs_Cone, GeomAbs_Sphere };
Standard_Integer isurf = 0;
if (n > 3)
isurf = Draw::Atoi(a[3]);
isurf = Min(isurf, 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]);
@@ -1079,6 +1096,7 @@ static Standard_Integer getanasurf(Draw_Interpretor& di,
if (!aRes.IsNull())
{
DrawTrSurf::Set(a[1], aRes);
di << "Gap = " << aCanonRec.GetGap() << "\n";
}
else
{
@@ -1094,6 +1112,12 @@ static Standard_Integer getanasurf(Draw_Interpretor& di,
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();
@@ -1106,11 +1130,18 @@ Standard_Integer getanacurve(Draw_Interpretor& di,
GeomAbs_CurveType aTargets[] = { GeomAbs_Line, GeomAbs_Circle, GeomAbs_Ellipse };
Standard_Integer icurv = 0;
if (n > 3)
icurv = Draw::Atoi(a[3]);
icurv = Min(icurv, 2);
Standard_Real tol = 1.e-7;
{
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]);
tol = Draw::Atof(a[4]);
ShapeAnalysis_CanonicalRecognition aCanonRec(sh);
Handle(Geom_Curve) aRes;
@@ -1144,6 +1175,7 @@ Standard_Integer getanacurve(Draw_Interpretor& di,
if (!aRes.IsNull())
{
DrawTrSurf::Set(a[1], aRes);
di << "Gap = " << aCanonRec.GetGap() << "\n";
}
else
{
@@ -1194,7 +1226,7 @@ Standard_Integer getanacurve(Draw_Interpretor& di,
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("getanasurf", "getanasurf res shape [target [tol [sample]]] ", __FILE__, getanasurf, g);
theCommands.Add("getanacurve", "getanacurve res shape [target [tol]]", __FILE__, getanacurve, g);

View File

@@ -47,3 +47,9 @@ 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

View File

@@ -46,6 +46,20 @@
#include <GeomAbs_SurfaceType.hxx>
#include <GeomAbs_CurveType.hxx>
#include <NCollection_Vector.hxx>
#include <GeomAdaptor_Curve.hxx>
#include <GeomAdaptor_Surface.hxx>
#include <BRepLib_FindSurface.hxx>
#include <TColgp_HArray1OfXYZ.hxx>
#include <math_Vector.hxx>
#include <ShapeAnalysis_FuncSphereLSDist.hxx>
#include <ShapeAnalysis_FuncCylinderLSDist.hxx>
#include <ShapeAnalysis_FuncConeLSDist.hxx>
#include <BRepAdaptor_Curve.hxx>
#include <GCPnts_AbscissaPoint.hxx>
#include <GCPnts_QuasiUniformAbscissa.hxx>
#include <math_PSO.hxx>
#include <math_Powell.hxx>
#include <ElSLib.hxx>
//Declaration of static functions
@@ -64,6 +78,17 @@ static Standard_Boolean CompareSurfParams(const GeomAbs_SurfaceType theTarget, c
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);
//=======================================================================
@@ -201,11 +226,27 @@ Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsPlane(const Standard_Real
TColStd_Array1OfReal aParams(1, 1);
//
GeomAbs_SurfaceType aTarget = GeomAbs_Plane;
Standard_Boolean isOK = IsElementarySurf(aTarget, theTol, aPos, aParams);
if (!isOK)
return isOK;
thePln.SetPosition(aPos);
return Standard_True;
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;
}
//=======================================================================
@@ -220,12 +261,41 @@ Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsCylinder(const Standard_R
aParams(1) = theCyl.Radius();
//
GeomAbs_SurfaceType aTarget = GeomAbs_Cylinder;
Standard_Boolean isOK = IsElementarySurf(aTarget, theTol, aPos, aParams);
if (!isOK)
return isOK;
theCyl.SetPosition(aPos);
theCyl.SetRadius(aParams(1));
return Standard_True;
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;
}
//=======================================================================
@@ -241,13 +311,44 @@ Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsCone(const Standard_Real
aParams(2) = theCone.RefRadius();
//
GeomAbs_SurfaceType aTarget = GeomAbs_Cone;
Standard_Boolean isOK = IsElementarySurf(aTarget, theTol, aPos, aParams);
if (!isOK)
return isOK;
theCone.SetPosition(aPos);
theCone.SetSemiAngle(aParams(1));
theCone.SetRadius(aParams(2));
return Standard_True;
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;
}
//=======================================================================
@@ -262,12 +363,41 @@ Standard_Boolean ShapeAnalysis_CanonicalRecognition::IsSphere(const Standard_Rea
aParams(1) = theSphere.Radius();
//
GeomAbs_SurfaceType aTarget = GeomAbs_Sphere;
Standard_Boolean isOK = IsElementarySurf(aTarget, theTol, aPos, aParams);
if (!isOK)
return isOK;
theSphere.SetPosition(aPos);
theSphere.SetRadius(aParams(1));
return Standard_True;
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;
}
//=======================================================================
@@ -557,14 +687,21 @@ Handle(Geom_Surface) ShapeAnalysis_CanonicalRecognition::GetSurface(const TopoDS
Standard_Integer ifit = -1, i;
Standard_Real aMinDev = RealLast();
for (i = 0; i < aSurfs.Size(); ++i)
if (aSurfs.Size() == 1)
{
SetSurfParams(theTarget, aSurfs(i), aPos, aParams);
Standard_Real aDev = DeviationSurfParams(theTarget, thePos, theParams, aPos, aParams);
if (aDev < aMinDev)
ifit = 0;
}
else
{
for (i = 0; i < aSurfs.Size(); ++i)
{
aMinDev = aDev;
ifit = 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)
@@ -649,6 +786,130 @@ Handle(Geom_Surface) ShapeAnalysis_CanonicalRecognition::GetSurface(const TopoDS
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;
}
//=======================================================================
@@ -749,37 +1010,32 @@ Standard_Integer GetNbPars(const GeomAbs_SurfaceType theTarget)
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)
{
Handle(Geom_Line) aConic = Handle(Geom_Line)::DownCast(theConic);
if (aConic.IsNull())
return Standard_False;
gp_Lin aLin = aConic->Lin();
gp_Lin aLin = aGAC.Line();
thePos.SetAxis(aLin.Position());
}
else if (theTarget == GeomAbs_Circle)
{
Handle(Geom_Circle) aConic = Handle(Geom_Circle)::DownCast(theConic);
if (aConic.IsNull())
return Standard_False;
gp_Circ aCirc = aConic->Circ();
gp_Circ aCirc = aGAC.Circle();
thePos = aCirc.Position();
theParams(1) = aCirc.Radius();
}
else if (theTarget == GeomAbs_Ellipse)
{
Handle(Geom_Ellipse) aConic = Handle(Geom_Ellipse)::DownCast(theConic);
if (aConic.IsNull())
return Standard_False;
gp_Elips anElips = aConic->Elips();
gp_Elips anElips = aGAC.Ellipse();
thePos = anElips.Position();
theParams(1) = anElips.MajorRadius();
theParams(2) = anElips.MinorRadius();
}
else
{
return Standard_False;
}
return Standard_True;
}
@@ -826,43 +1082,37 @@ Standard_Boolean SetSurfParams(const GeomAbs_SurfaceType theTarget, const Handle
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)
{
Handle(Geom_Plane) aSurf = Handle(Geom_Plane)::DownCast(theElemSurf);
if (aSurf.IsNull())
return Standard_False;
gp_Pln aPln = aSurf->Pln();
gp_Pln aPln = aGAS.Plane();
thePos = aPln.Position();
}
else if (theTarget == GeomAbs_Cylinder)
{
Handle(Geom_CylindricalSurface) aSurf = Handle(Geom_CylindricalSurface)::DownCast(theElemSurf);
if (aSurf.IsNull())
return Standard_False;
gp_Cylinder aCyl = aSurf->Cylinder();
gp_Cylinder aCyl = aGAS.Cylinder();
thePos = aCyl.Position();
theParams(1) = aCyl.Radius();
}
else if (theTarget == GeomAbs_Cone)
{
Handle(Geom_ConicalSurface) aSurf = Handle(Geom_ConicalSurface)::DownCast(theElemSurf);
if (aSurf.IsNull())
return Standard_False;
gp_Cone aCon = aSurf->Cone();
gp_Cone aCon = aGAS.Cone();
thePos = aCon.Position();
theParams(1) = aCon.SemiAngle();
theParams(2) = aCon.RefRadius();
}
else if (theTarget == GeomAbs_Sphere)
{
Handle(Geom_SphericalSurface) aSurf = Handle(Geom_SphericalSurface)::DownCast(theElemSurf);
if (aSurf.IsNull())
return Standard_False;
gp_Sphere aSph = aSurf->Sphere();
gp_Sphere aSph = aGAS.Sphere();
thePos = aSph.Position();
theParams(1) = aSph.Radius();
}
@@ -962,3 +1212,204 @@ Standard_Real DeviationSurfParams(const GeomAbs_SurfaceType theTarget,
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<Standard_Real> aLengths;
NCollection_Vector<BRepAdaptor_Curve> aCurves;
NCollection_Vector<gp_XYZ> 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
}
}

View File

@@ -150,6 +150,11 @@ private:
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:

View File

@@ -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 <ShapeAnalysis_FuncConeLSDist.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include <gp_Ax3.hxx>
#include <math_Vector.hxx>
#include <ElSLib.hxx>
//=======================================================================
//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;
}

View File

@@ -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 <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
#include <math_MultipleVarFunction.hxx>
#include <TColgp_HArray1OfXYZ.hxx>
#include <math_Vector.hxx>
#include <gp_Dir.hxx>
//! 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

View File

@@ -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 <ShapeAnalysis_FuncCylinderLSDist.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include <math_Vector.hxx>
//=======================================================================
//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;
}

View File

@@ -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 <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
#include <math_MultipleVarFunctionWithGradient.hxx>
#include <TColgp_HArray1OfXYZ.hxx>
#include <math_Vector.hxx>
#include <gp_Dir.hxx>
//! 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

View File

@@ -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 <ShapeAnalysis_FuncSphereLSDist.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include <math_Vector.hxx>
//=======================================================================
//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;
}

View File

@@ -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 <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
#include <math_MultipleVarFunctionWithGradient.hxx>
#include <TColgp_HArray1OfXYZ.hxx>
#include <math_Vector.hxx>
//! 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

View File

@@ -54,7 +54,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);
@@ -77,7 +77,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);
@@ -86,14 +86,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()) {
@@ -140,7 +147,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;
}
}

12
tests/cr/approx/A1 Normal file
View File

@@ -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"
}

15
tests/cr/approx/A2 Normal file
View File

@@ -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"
}

15
tests/cr/approx/A3 Normal file
View File

@@ -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"
}

16
tests/cr/approx/A4 Normal file
View File

@@ -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"
}

View File

@@ -1,6 +1,6 @@
beziercurve bz 3 0 0 0 .5 .5e-3 0 1 0 0
mkedge e bz
getanacurve c e 0 1.e-3
getanacurve c e lin 1.e-3
if {[isdraw c]} {
set log [dump c]
if { [regexp {Line} $log ] != 1 } {

View File

@@ -1,6 +1,6 @@
beziercurve bz 3 0 0 0 .5 .5e-3 0 1 0 0
mkedge e bz
getanacurve c e 1 1.e-3
getanacurve c e cir 1.e-3
if {[isdraw c]} {
set log [dump c]
if { [regexp {Circle} $log ] != 1 } {

View File

@@ -1,7 +1,7 @@
ellipse bz 0 0 0 0 0 1 1 .5
convert bz bz
mkedge e bz
getanacurve c e 2 1.e-7
getanacurve c e ell 1.e-7
if {[isdraw c]} {
set log [dump c]
if { [regexp {Ellipse} $log ] != 1 } {

View File

@@ -3,7 +3,7 @@ mkedge e1 bz 0 .3
mkedge e2 bz .3 0.7
mkedge e3 bz .7 1.
wire w e1 e2 e3
getanacurve c w 0 1.e-3
getanacurve c w lin 1.e-3
if {[isdraw c]} {
set log [dump c]
if { [regexp {Line} $log ] != 1 } {

View File

@@ -4,7 +4,7 @@ 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 1 1.e-7
getanacurve c w cir 1.e-7
if {[isdraw c]} {
set log [dump c]
if { [regexp {Circle} $log ] != 1 } {

View File

@@ -4,7 +4,7 @@ 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 2 1.e-7
getanacurve c w ell 1.e-7
if {[isdraw c]} {
set log [dump c]
if { [regexp {Ellipse} $log ] != 1 } {

View File

@@ -2,7 +2,7 @@ plane surf
trim surf surf -1 1 -1 1
convert surf surf
mkface f surf
getanasurf asurf f 0 1.e-7
getanasurf asurf f pln 1.e-7
if {[isdraw asurf]} {
set log [dump asurf]
if { [regexp {Plane} $log ] != 1 } {

View File

@@ -2,7 +2,7 @@ cylinder surf 1
trimv surf surf -1 1
convert surf surf
mkface f surf
getanasurf asurf f 1 1.e-7
getanasurf asurf f cyl 1.e-7
if {[isdraw asurf]} {
set log [dump asurf]
if { [regexp {CylindricalSurface} $log ] != 1 } {

View File

@@ -2,7 +2,7 @@ cone surf 30 0
trimv surf surf -1 0
convert surf surf
mkface f surf
getanasurf asurf f 2 1.e-7
getanasurf asurf f con 1.e-7
if {[isdraw asurf]} {
set log [dump asurf]
if { [regexp {ConicalSurface} $log ] != 1 } {

View File

@@ -1,7 +1,7 @@
sphere surf 1
convert surf surf
mkface f surf
getanasurf asurf f 3 1.e-7
getanasurf asurf f sph 1.e-7
if {[isdraw asurf]} {
set log [dump asurf]
if { [regexp {SphericalSurface} $log ] != 1 } {

View File

@@ -12,7 +12,7 @@ trim surf4 surf 0 1 -1 0
convert surf4 surf4
mkface f4 surf4
sewing sh f1 f2 f3 f4
getanasurf asurf sh 0 1.e-7
getanasurf asurf sh pln 1.e-7
if {[isdraw asurf]} {
set log [dump asurf]
if { [regexp {Plane} $log ] != 1 } {

View File

@@ -12,7 +12,7 @@ trim surf4 surf 3 6 -1 0
convert surf4 surf4
mkface f4 surf4
sewing sh f1 f2 f3 f4
getanasurf asurf sh 1 1.e-7
getanasurf asurf sh cyl 1.e-7
if {[isdraw asurf]} {
set log [dump asurf]
if { [regexp {CylindricalSurface} $log ] != 1 } {

View File

@@ -12,7 +12,7 @@ trim surf4 surf 3 6 0 1
convert surf4 surf4
mkface f4 surf4
sewing sh f1 f2 f3 f4
getanasurf asurf sh 2 1.e-7
getanasurf asurf sh con 1.e-7
if {[isdraw asurf]} {
set log [dump asurf]
if { [regexp {ConicalSurface} $log ] != 1 } {

View File

@@ -12,7 +12,7 @@ trim surf4 surf 3 6 -1.5 0
convert surf4 surf4
mkface f4 surf4
sewing sh f1 f2 f3 f4
getanasurf asurf sh 3 1.e-7
getanasurf asurf sh sph 1.e-7
if {[isdraw asurf]} {
set log [dump asurf]
if { [regexp {SphericalSurface} $log ] != 1 } {

View File

@@ -24,7 +24,7 @@ add ss_2 w
add ss_3 w
add ss_4 w
add ss_5 w
getanasurf asurf w 1 1.e-7 cyl
getanasurf asurf w cyl 1.e-7 cyl
if {[isdraw asurf]} {
set log [dump asurf]
if { [regexp {CylindricalSurface} $log ] != 1 } {

View File

@@ -1 +1,2 @@
001 base
002 approx