--- /dev/null
+# Curvature of a Face along given direction
+
+import salome
+salome.salome_init_without_session()
+import GEOM
+from salome.geom import geomBuilder
+geompy = geomBuilder.New()
+import math
+
+O = geompy.MakeVertex(0, 0, 0, 'O')
+OX = geompy.MakeVectorDXDYDZ(1, 0, 0, 'OX')
+OY = geompy.MakeVectorDXDYDZ(0, 1, 0, 'OY')
+OZ = geompy.MakeVectorDXDYDZ(0, 0, 1, 'OZ')
+
+pXYZ = geompy.MakeVertex(105, 105, 105, 'pXYZ')
+pY = geompy.MakeVertex(0, 105, 0, 'pY')
+pZ = geompy.MakeVertex(0, 0, 105, 'pZ')
+
+vZ_XY = geompy.MakeVectorDXDYDZ(-1, -1, 1, 'vZ-XY')
+vZ_XY2 = geompy.MakeVectorDXDYDZ(-1, -1, 10, 'vZ-XY')
+vZ_XY3 = geompy.MakeVectorDXDYDZ(-1, -1, 100, 'vZ-XY')
+
+R = 100.0
+
+# I. Curvature of a Sphere
+Sphere_1 = geompy.MakeSphereR(R, 'Sphere_1')
+[Sph] = geompy.ExtractShapes(Sphere_1, geompy.ShapeType["FACE"], True, "Sph")
+
+curvature_1 = geompy.CurvatureOnFace(Sph, pXYZ, OX, 'curvature_sph_pXYZ_OX')
+curvature_2 = geompy.CurvatureOnFace(Sph, pXYZ, vZ_XY, 'curvature_sph_pXYZ_vt')
+curvature_3 = geompy.CurvatureOnFace(Sph, pY, OX, 'curvature_sph_pY_OX')
+
+# All sphere curvature radiuces = R
+assert(abs(geompy.BasicProperties(curvature_1)[0] - R) < 1e-07)
+assert(abs(geompy.BasicProperties(curvature_2)[0] - R) < 1e-07)
+assert(abs(geompy.BasicProperties(curvature_3)[0] - R) < 1e-07)
+
+# Normal direction
+isExcept = False
+try:
+ geompy.CurvatureOnFace(Sph, pY, OY)
+except:
+ isExcept = True
+assert(isExcept)
+
+# Pole (min and max curvatures are not defined, find via line projection?)
+isExcept = False
+try:
+ geompy.CurvatureOnFace(Sph, pZ, OX, 'curvature_sph_pZ_OX')
+except:
+ isExcept = True
+ print(geompy.MeasuOp.GetErrorCode())
+assert(isExcept)
+
+# II. Curvature of a Cylinder
+Cylinder_1 = geompy.MakeCylinderRH(R, 300, 'Cylinder_1')
+[Face_1,Face_2,Face_3] = geompy.ExtractShapes(Cylinder_1, geompy.ShapeType["FACE"], True, "Face")
+
+# Curvature radius of a cylinder along any direction, orthogonal to its Z axis, equal to R
+curvature_4 = geompy.CurvatureOnFace(Face_2, pY, OX, 'curvature_cyl_pY_OX')
+assert(abs(geompy.BasicProperties(curvature_4)[0] - R) < 1e-07)
+
+# Curvature radius of a cylinder along its Z direction is infinite
+curvature_zero = geompy.CurvatureOnFace(Face_2, pY, OZ)
+assert(geompy.MeasuOp.GetErrorCode() == "ZERO_CURVATURE")
+assert(not curvature_zero)
+
+# Curvature radius of a cylinder along some direction, different from two above
+curvature_5 = geompy.CurvatureOnFace(Face_2, pY, vZ_XY, 'curvature_cyl_pY_vZ_XY')
+curvature_6 = geompy.CurvatureOnFace(Face_2, pY, vZ_XY2, 'curvature_cyl_pY_vZ_XY2')
+curvature_7 = geompy.CurvatureOnFace(Face_2, pY, vZ_XY3, 'curvature_cyl_pY_vZ_XY3')
+
+# R < r5 < r6 < r7
+# r5 = 100.01, r6 = 101.0, r7 = 200
+r5 = geompy.BasicProperties(curvature_5)[0]
+r6 = geompy.BasicProperties(curvature_6)[0]
+r7 = geompy.BasicProperties(curvature_7)[0]
+
+assert(R + 1e-07 < r5)
+assert(r5 + 1e-07 < r6)
+assert(r6 + 1e-07 < r7)
+
+# Projection aborted. Point is out of the face boundaries.
+isExcept = False
+try:
+ pXY_Z = geompy.MakeVertex(105, 105, -105, 'pXY_Z')
+ geompy.CurvatureOnFace(Face_2, pXY_Z, OX, 'curvature_cyl_pXY_Z')
+except:
+ isExcept = True
+assert(isExcept)
+
+# Projection aborted (point on axis). Equal distances to many points.
+isExcept = False
+try:
+ geompy.CurvatureOnFace(Face_2, O, vZ_XY, 'curvature_cyl_O')
+except:
+ isExcept = True
+assert(isExcept)
+
+# Curvature radius of a planar face is infinite
+curvature_zero_2 = geompy.CurvatureOnFace(Face_1, pZ, OX)
+assert(geompy.MeasuOp.GetErrorCode() == "ZERO_CURVATURE")
+assert(not curvature_zero_2)
+
+# III. Curvature of a "Horse saddle"
+[Edge_1,Edge_2,Edge_3] = geompy.ExtractShapes(Sphere_1, geompy.ShapeType["EDGE"], True)
+geompy.addToStudyInFather( Sphere_1, Edge_1, 'Edge_1' )
+geompy.addToStudyInFather( Sphere_1, Edge_2, 'Edge_2' )
+geompy.addToStudyInFather( Sphere_1, Edge_3, 'Edge_3' )
+
+Rotation_1 = geompy.MakeRotation(Edge_3, OX, 90*math.pi/180.0, 'Rotation_1')
+Rotation_2 = geompy.MakeRotation(Rotation_1, OY, 180*math.pi/180.0, 'Rotation_2')
+Translation_1 = geompy.MakeTranslation(Rotation_2, 200, 0, 0, 'Translation_1')
+Translation_2 = geompy.MakeTranslation(Edge_3, 100, 100, 0, 'Translation_2')
+Translation_3 = geompy.MakeTranslation(Translation_2, 0, -200, 0, 'Translation_3')
+Filling_1 = geompy.MakeFilling([Translation_2, Edge_3, Translation_3])
+geompy.addToStudy(Filling_1, 'Filling_1')
+Vertex_2 = geompy.MakeVertex(100, 0, 0, 'Vertex_2')
+
+curvature_Y = geompy.CurvatureOnFace(Filling_1, Vertex_2, OY, 'curvature_Y')
+curvature_Z = geompy.CurvatureOnFace(Filling_1, Vertex_2, OZ, 'curvature_Z')
+
+cury = geompy.VectorCoordinates(curvature_Y)
+curz = geompy.VectorCoordinates(curvature_Z)
+
+# Vectors should be opposite, scalar product should be negative
+assert(cury[0]*curz[0] + cury[1]*curz[1] + cury[2]*curz[2] < -1e-07)
+
+# Normal direction
+norm_1 = geompy.GetNormal(Filling_1, Vertex_2, "Normal_1")
+isExcept = False
+try:
+ geompy.CurvatureOnFace(Filling_1, Vertex_2, norm_1)
+except:
+ isExcept = True
+assert(isExcept)
import_export.py
inertia.py
min_distance.py
+ curvature_face.py
normal_face.py
notebook_geom.py
polyline.py
\until geompy.GetSubShapesWithTolerance(Box, GEOM.FACE, GEOM.CC_LE, 1.e-7, "le")
\anchor swig_MakeExtraction
-\until print "DONE"
+\until geompy.MakeExtraction(Box, [16], "Ext_no_vertex")
+
+\anchor swig_CurvatureOnFace
+\until print("DONE")
*/
<source>STB_NORMALE</source>
<translation>Compute normal to the face</translation>
</message>
+ <message>
+ <source>MEN_CURVATURE</source>
+ <translation>Vector of curvature</translation>
+ </message>
<message>
<source>TOP_MEASURE_ANGLE</source>
<translation>Angle</translation>
<source>STB_NORMALE</source>
<translation>Vecteur normal à une face</translation>
</message>
+ <message>
+ <source>MEN_CURVATURE</source>
+ <translation>Vecteur de courbure</translation>
+ </message>
<message>
<source>TOP_MEASURE_ANGLE</source>
<translation>Angle</translation>
<source>STB_NORMALE</source>
<translation>フェースに垂直</translation>
</message>
+ <message>
+ <source>MEN_CURVATURE</source>
+ <translation>Vector_of_curvature</translation>
+ </message>
<message>
<source>TOP_MEASURE_ANGLE</source>
<translation>角度</translation>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeVertex.hxx>
#include <BRepPrimAPI_MakeBox.hxx>
+#include <BRepOffsetAPI_NormalProjection.hxx>
#include <TopAbs.hxx>
#include <TopoDS.hxx>
#include <GProp_GProps.hxx>
#include <GeomLProp_SLProps.hxx>
+#include <GeomLProp_CLProps.hxx>
#include <Geom_Surface.hxx>
+#include <GeomAPI_ProjectPointOnCurve.hxx>
#include <Bnd_Box.hxx>
#include <BRepBndLib.hxx>
const TopoDS_Shape& thePoint,
const TopoDS_Shape& theDir)
{
- if (theFace.IsNull())
- Standard_NullObject::Raise("Face for curvature calculation is null");
- if (theFace.ShapeType() != TopAbs_FACE)
- Standard_NullObject::Raise("Shape for curvature calculation is not a face");
- TopoDS_Face aFace = TopoDS::Face(theFace);
- Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
- if (aSurf.IsNull())
- Standard_NullObject::Raise("Surface for curvature calculation is null");
-
+ // Point
if (thePoint.IsNull())
Standard_NullObject::Raise("Point for curvature measurement is null");
if (thePoint.ShapeType() != TopAbs_VERTEX)
- Standard_NullObject::Raise("Point for curvature calculation is not a vertex");
+ Standard_TypeMismatch::Raise("Point for curvature calculation is not a vertex");
gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(thePoint));
- gp_Vec aV = GEOMUtils::GetVector(theDir, Standard_False);
+ // Point projection on the face
+ Standard_Real U, V;
+ aPnt = GEOMUtils::ProjectPointOnFace(aPnt, theFace, U, V);
+ gp_Pnt2d UV (U, V);
- // Point projection parameters on surface
- ShapeAnalysis_Surface aSAS (aSurf);
- gp_Pnt2d UV = aSAS.ValueOfUV(aPnt, Precision::Confusion());
- aPnt = aSAS.Value(UV);
+ // Face and Vector
+ TopoDS_Face aFace = TopoDS::Face(theFace);
+ gp_Vec aV = GEOMUtils::GetVector(theDir, Standard_False);
// Calculate differential properties
BRepAdaptor_Surface aSurfAdapt (aFace);
BRepLProp_SLProps Props (aSurfAdapt, UV.X(), UV.Y(), 2, 1e-7);
if (!Props.IsCurvatureDefined())
- Standard_NullObject::Raise("Curvature calculation failed");
+ Standard_ConstructionError::Raise("Curvature calculation failed");
// Get differential properties
gp_Vec Xu = Props.D1U();
gp_Vec2d T (aV.Dot(aDirU), aV.Dot(aDirV));
if (Abs(T.X()) < Precision::Confusion() &&
Abs(T.Y()) < Precision::Confusion())
- Standard_NullObject::Raise("Curvature calculation failed: direction is normal to the face");
+ Standard_ConstructionError::Raise
+ ("Curvature calculation failed: direction is normal to the face");
// Coefficients of the FFF
double E = Xu.Dot(Xu);
double M = n.Dot(Xuv);
double N = n.Dot(Xvv);
- // Calculate curvature using the coefficients of both fundamental forms
- double k = 0.;
+ // Calculate radius (or -radius) of curvature
+ // using the coefficients of both fundamental forms
+ double r = 0.;
if (Abs(T.X()) < Precision::Confusion()) {
- //if (Abs(G) < Precision::Confusion())
- // Standard_NullObject::Raise("Curvature calculation failed: G = 0");
- //k = N / G; // curvature
if (Abs(N) < Precision::Confusion())
- Standard_NullObject::Raise("Curvature calculation failed: N = 0");
- k = G / N; // radius of curvature
+ Standard_Failure::Raise("ZERO_CURVATURE");
+ r = G / N;
}
else {
double lambda = T.Y() / T.X();
double detE = E + 2*F*lambda + G*lambda*lambda;
double detL = L + 2*M*lambda + N*lambda*lambda;
- //if (Abs(detE) < Precision::Confusion())
if (Abs(detL) < Precision::Confusion())
- Standard_NullObject::Raise("Curvature calculation failed: det = 0");
- //k = detL / detE; // curvature
- k = detE / detL; // radius of curvature
+ Standard_Failure::Raise("ZERO_CURVATURE");
+ r = detE / detL;
}
// Result
gp_Dir aNormal (n);
- gp_Pnt aPntEnd (aPnt.XYZ() + aNormal.XYZ() * k);
+ gp_Pnt aPntEnd (aPnt.XYZ() + aNormal.XYZ() * r);
BRepBuilderAPI_MakeEdge aBuilder (aPnt, aPntEnd);
if (!aBuilder.IsDone())
- Standard_NullObject::Raise("Curvature calculation failed: edge is not built");
+ Standard_ConstructionError::Raise("Curvature calculation failed: edge is not built");
return aBuilder.Shape();
}
AddParam( theParams, "Face", aCI.GetBase() );
AddParam( theParams, "Point", aCI.GetPoint(), "face center" );
break;
+ case CURVATURE_VEC_MEASURE:
+ theOperationName = "CURVATURE";
+ AddParam( theParams, "Face", aCI.GetBase() );
+ AddParam( theParams, "Point", aCI.GetPoint(), "point of interest" );
+ AddParam( theParams, "Vector", aCI.GetDirection(), "direction of interest" );
+ break;
default:
return false;
}
Handle(GEOM_Function) aTargetFunction = TI.GetPlane();
if (aTargetFunction.IsNull()) return 0;
TopoDS_Shape aFaceShape = aTargetFunction->GetValue();
- //if (aFaceShape.IsNull() || aFaceShape.ShapeType() != TopAbs_FACE) {
- // Standard_ConstructionError::Raise
- // ("Projection aborted : the target shape is not a face");
- //}
-
- Standard_Real tol = 1.e-4;
if (anOriginal.ShapeType() == TopAbs_VERTEX) {
- if (aFaceShape.IsNull() || aFaceShape.ShapeType() != TopAbs_FACE) {
- Standard_ConstructionError::Raise
- ("Projection aborted : the target shape is not a face");
- }
- TopoDS_Face aFace = TopoDS::Face(aFaceShape);
- Handle(Geom_Surface) surface = BRep_Tool::Surface(aFace);
- double U1, U2, V1, V2;
- //surface->Bounds(U1, U2, V1, V2);
- BRepTools::UVBounds(aFace, U1, U2, V1, V2);
-
- // projector
- GeomAPI_ProjectPointOnSurf proj;
- proj.Init(surface, U1, U2, V1, V2, tol);
-
- gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(anOriginal));
- proj.Perform(aPnt);
- if (!proj.IsDone()) {
- Standard_ConstructionError::Raise
- ("Projection aborted : the algorithm failed");
- }
- int nbPoints = proj.NbPoints();
- if (nbPoints < 1) {
- Standard_ConstructionError::Raise("No solution found");
- }
-
Standard_Real U, V;
- proj.LowerDistanceParameters(U, V);
- gp_Pnt2d aProjPnt (U, V);
-
- // classifier
- BRepClass_FaceClassifier aClsf (aFace, aProjPnt, tol);
- if (aClsf.State() != TopAbs_IN && aClsf.State() != TopAbs_ON) {
- bool isSol = false;
- double minDist = RealLast();
- for (int i = 1; i <= nbPoints; i++) {
- Standard_Real Ui, Vi;
- proj.Parameters(i, Ui, Vi);
- aProjPnt = gp_Pnt2d(Ui, Vi);
- aClsf.Perform(aFace, aProjPnt, tol);
- if (aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON) {
- isSol = true;
- double dist = proj.Distance(i);
- if (dist < minDist) {
- minDist = dist;
- U = Ui;
- V = Vi;
- }
- }
- }
- if (!isSol) {
- Standard_ConstructionError::Raise("No solution found");
- }
- }
-
- gp_Pnt surfPnt = surface->Value(U, V);
-
+ gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(anOriginal));
+ gp_Pnt surfPnt = GEOMUtils::ProjectPointOnFace(aPnt, aFaceShape, U, V);
aShape = BRepBuilderAPI_MakeVertex(surfPnt).Shape();
}
else {
#include <BRepGProp.hxx>
#include <BRepTools.hxx>
+#include <BRepClass_FaceClassifier.hxx>
#include <BRepClass3d_SolidClassifier.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <TopTools_ListIteratorOfListOfShape.hxx>
#include <TopTools_Array1OfShape.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
+
#include <Geom_Circle.hxx>
#include <Geom_Surface.hxx>
#include <Geom_Plane.hxx>
return aResult;
}
+//=======================================================================
+// function : ProjectPointOnFace()
+// purpose : Returns the projection (3d point) if found, throws an exception otherwise
+//=======================================================================
+gp_Pnt GEOMUtils::ProjectPointOnFace(const gp_Pnt& thePoint,
+ const TopoDS_Shape& theFace,
+ double& theU, double& theV)
+{
+ if (theFace.IsNull() || theFace.ShapeType() != TopAbs_FACE)
+ Standard_TypeMismatch::Raise
+ ("Projection aborted : the target shape is not a face");
+
+ TopoDS_Face aFace = TopoDS::Face(theFace);
+ Handle(Geom_Surface) surface = BRep_Tool::Surface(aFace);
+ double U1, U2, V1, V2;
+ BRepTools::UVBounds(aFace, U1, U2, V1, V2);
+
+ // projector
+ Standard_Real tol = 1.e-4;
+ GeomAPI_ProjectPointOnSurf proj;
+ proj.Init(surface, U1, U2, V1, V2, tol);
+ proj.Perform(thePoint);
+ if (!proj.IsDone())
+ StdFail_NotDone::Raise("Projection aborted : the algorithm failed");
+ int nbPoints = proj.NbPoints();
+ if (nbPoints < 1)
+ Standard_ConstructionError::Raise("Projection aborted : No solution found");
+ proj.LowerDistanceParameters(theU, theV);
+ gp_Pnt2d aProjPnt (theU, theV);
+
+ // classifier
+ BRepClass_FaceClassifier aClsf (aFace, aProjPnt, tol);
+ if (aClsf.State() != TopAbs_IN && aClsf.State() != TopAbs_ON) {
+ bool isSol = false;
+ double minDist = RealLast();
+ for (int i = 1; i <= nbPoints; i++) {
+ Standard_Real Ui, Vi;
+ proj.Parameters(i, Ui, Vi);
+ aProjPnt = gp_Pnt2d(Ui, Vi);
+ aClsf.Perform(aFace, aProjPnt, tol);
+ if (aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON) {
+ isSol = true;
+ double dist = proj.Distance(i);
+ if (dist < minDist) {
+ minDist = dist;
+ theU = Ui;
+ theV = Vi;
+ }
+ }
+ }
+ if (!isSol) {
+ Standard_ConstructionError::Raise("Projection aborted : No solution found");
+ }
+ }
+
+ gp_Pnt surfPnt = surface->Value(theU, theV);
+ return surfPnt;
+}
+
//=======================================================================
// function : ConvertClickToPoint()
// purpose : Returns the point clicked in 3D view
const TopoDS_Shape& theShape2,
gp_Pnt& thePnt1, gp_Pnt& thePnt2);
+ /*!
+ * \brief Computes normal projection of \a thePoint to \a theFace.
+ *
+ * \param thePoint the 3d point
+ * \param theFace the face shape
+ * \param theU the output U parameter of the point on the face
+ * \param theV the output V parameter of the point on the face
+ * \retval the projection (3d point) if found, throws an exception otherwise
+ */
+ Standard_EXPORT gp_Pnt ProjectPointOnFace(const gp_Pnt& thePoint,
+ const TopoDS_Shape& theFace,
+ double& theU, double& theV);
+
/*!
* \brief Returns the point clicked in 3D view.
*
geompy.MakeExtraction(Box, [18], "Ext_no_edge")
geompy.MakeExtraction(Box, [16], "Ext_no_vertex")
+ # CurvatureOnFace
+ Cylinder_1 = geompy.MakeCylinderRH(100, 50, 'Cylinder_r100_h150')
+ [Face_1,Face_2,Face_3] = geompy.ExtractShapes(Cylinder_1, geompy.ShapeType["FACE"], True, "Face")
+ curvature_1 = geompy.CurvatureOnFace(Face_2, px, vy, 'curvature_cyl_px_vy')
+ assert(abs(geompy.BasicProperties(curvature_1)[0] - 100) < 1e-07)
+ isExcept = False
+ try:
+ geompy.CurvatureOnFace(Face_2, px, vz)
+ except:
+ isExcept = True
+ assert(geompy.MeasuOp.GetErrorCode() == "Curvature radius is infinite")
+ assert(isExcept)
print("DONE")
return aSurf
## @}
- ## Measure curvature of surface in the given point along the given direction.
+ ## Measure curvature radius of surface in the given point along the given direction.
# @param theSurf the given face.
# @param thePoint given point.
# @param theDirection given direction.
# in the given point in case of positive curvature value
# and opposite to the normal in case of negative curvature.
# The normal of the returned vector is equal to the
- # absolute value of the curvature.
+ # absolute value of the curvature radius.
+ # Null shape is returned in case of infinite radius
+ # (zero curvature), for example, in case of flat face.
#
- ## @ref swig_todo "Example"
+ ## @ref swig_CurvatureOnFace "Example"
@ManageTransactions("MeasuOp")
def CurvatureOnFace(self, theSurf, thePoint, theDirection, theName=None):
"""
- Measure curvature of surface in the given point along the given direction.
+ Measure curvature radius of surface in the given point along the given direction.
Parameters:
theSurf the given face.
in the given point in case of positive curvature value
and opposite to the normal in case of negative curvature.
The normal of the returned vector is equal to the
- absolute value of the curvature.
+ absolute value of the curvature radius.
+ Null shape is returned in case of infinite radius
+ (zero curvature), for example, in case of flat face.
Example of usage:
curvature_1 = geompy.CurvatureOnFace(Face_1, Vertex_1, OX)
"""
aVec = self.MeasuOp.SurfaceCurvatureByPointAndDirection(theSurf,thePoint,theDirection)
- RaiseIfFailed("CurvatureOnFace", self.MeasuOp)
- self._autoPublish(aVec, theName, "curvature")
+ if self.MeasuOp.GetErrorCode() != "ZERO_CURVATURE":
+ RaiseIfFailed("CurvatureOnFace", self.MeasuOp)
+ self._autoPublish(aVec, theName, "curvature")
return aVec
## Get min and max tolerances of sub-shapes of theShape
