[bos #29472] [EDF] (2022-T1) Advanced geometry features: curvature vector on a point of a face

doc/salome/examples/curvature_face.py

@@ -0,0 +1,203 @@
+# 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
+import numpy as np
+
+def test_acceptance():
+  """
+  Acceptance test [tuleap29472]
+  """
+  Vector = [0,100,100]
+  O = geompy.MakeVertex(0, 0, 0)
+  OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
+  OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
+  OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
+  Cylinder_1 = geompy.MakeCylinderRH(100, 300)
+  Translation_1 = geompy.MakeTranslation(Cylinder_1, 0, 0, -150)
+  Vertex_1 = geompy.MakeVertex(100, 0, 0)
+  Vertex_2 = geompy.MakeVertex(100, -Vector[2], Vector[1])
+  Line_1 = geompy.MakeLineTwoPnt(Vertex_1, Vertex_2)
+  Plane_1 = geompy.MakePlane(Vertex_1, Line_1, 2000)
+  Rotation_1 = geompy.MakeRotation(Translation_1, OZ, 90*math.pi/180.0)# avoid to have degenerated edge across Vertex_1
+
+  [Face_1,Face_2,Face_3] = geompy.ExtractShapes(Rotation_1, geompy.ShapeType["FACE"], True)
+
+  curvature_29472 = np.array( geompy.VectorCoordinates( geompy.CurvatureOnFace(Face_2, Vertex_1, geompy.MakeVectorDXDYDZ(*Vector))) ).reshape(1,3)
+  expected_curvature = np.array( [-200.0,0.0,0.0] ).reshape(1,3)
+  assert( np.isclose( 0.0, np.linalg.norm( curvature_29472 - expected_curvature ) ,rtol=0,atol=1e-5 ) )
+
+  Intersection_1 = geompy.MakeSection(Face_2, Plane_1, True)
+  geompy.addToStudy( O, 'O' )
+  geompy.addToStudy( OX, 'OX' )
+  geompy.addToStudy( OY, 'OY' )
+  geompy.addToStudy( OZ, 'OZ' )
+  geompy.addToStudy( Vertex_1, 'Vertex_1' )
+  geompy.addToStudy( Cylinder_1, 'Cylinder_1' )
+  geompy.addToStudy( Translation_1, 'Translation_1' )
+  geompy.addToStudy( Vertex_2, 'Vertex_2' )
+  geompy.addToStudy( Line_1, 'Line_1' )
+  geompy.addToStudy( Plane_1, 'Plane_1' )
+  geompy.addToStudy( Rotation_1, 'Rotation_1' )
+  geompy.addToStudyInFather( Rotation_1, Face_1, 'Face_1' )
+  geompy.addToStudyInFather( Rotation_1, Face_2, 'Face_2' )
+  geompy.addToStudyInFather( Rotation_1, Face_3, 'Face_3' )
+  geompy.addToStudy( Intersection_1, 'Intersection_1' )
+  angle = math.asin(Vector[2]/math.sqrt(Vector[1]*Vector[1]+Vector[2]*Vector[2]))
+  tmp = geompy.MakeTranslation(Intersection_1,*[-elt for elt in geompy.PointCoordinates(Vertex_1)])
+  tmp = geompy.MakeRotation(tmp,OX,-angle)
+  Intersection_1_OXY = geompy.MakeTranslation(tmp,*geompy.PointCoordinates(Vertex_1))
+  geompy.addToStudy( Intersection_1_OXY, 'Intersection_1_OXY' )
+
+  eps = 0.01
+  offset = 0.75
+  p0 = np.array( geompy.PointCoordinates( geompy.MakeVertexOnCurve(Intersection_1_OXY,offset-eps) ) ).reshape(1,3)
+  p1 = np.array( geompy.PointCoordinates( geompy.MakeVertexOnCurve(Intersection_1_OXY,offset) ) ).reshape(1,3)
+  p2 = np.array( geompy.PointCoordinates( geompy.MakeVertexOnCurve(Intersection_1_OXY,offset+eps) ) ).reshape(1,3)
+  assert( np.isclose(0.0,np.linalg.norm(p1- np.array(geompy.PointCoordinates(Vertex_1)).reshape(1,3)  ),rtol=0,atol=1e-8) )
+  p01=(p0+p1)/2
+  p12=(p1+p2)/2
+  v0 = (p1-p0)/np.linalg.norm(p1-p0)
+  v1 = (p2-p1)/np.linalg.norm(p2-p1)
+  computedRadius =  1/np.linalg.norm((v1-v0)/np.linalg.norm(p12-p01))
+  # manual detection of radius : https://fr.wikipedia.org/wiki/Courbure_d%27un_arc
+  circle = geompy.MakeCircle(O,OZ,computedRadius)
+  circle = geompy.MakeTranslation(circle,100-computedRadius,0,0)
+  geompy.addToStudy(circle, "expectedCircle")
+  print("Radius expected is {}".format(computedRadius))
+  print("Radius obtain by CurvatureOnFace is {}".format(np.linalg.norm(curvature_29472)))
+
+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)
+
+# Pole
+isExcept = False
+try:
+  geompy.CurvatureOnFace(Sph, pZ, OX)
+except:
+  isExcept = True
+assert(isExcept)
+
+# Normal direction
+isExcept = False
+try:
+  geompy.CurvatureOnFace(Sph, pY,   OY)
+except:
+  isExcept = True
+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 = np.array( geompy.VectorCoordinates(curvature_Y) ).reshape(1,3)
+curz = np.array( geompy.VectorCoordinates(curvature_Z) ).reshape(1,3)
+cury_expected = np.array( [50,0,0] ).reshape(1,3)
+curz_expected = np.array( [-100,0,0] ).reshape(1,3)
+assert( np.isclose( 0.0, np.linalg.norm( cury - cury_expected ) ,rtol=0,atol=1e-5 ) )
+assert( np.isclose( 0.0, np.linalg.norm( curz - curz_expected ) ,rtol=0,atol=1e-5 ) )
+
+# 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)
+
+# acceptance case
+test_acceptance()

doc/salome/examples/tests.set

@@ -77,6 +77,7 @@ SET(GOOD_TESTS
   import_export.py
   inertia.py
   min_distance.py
+  curvature_face.py
   normal_face.py
   notebook_geom.py
   polyline.py

doc/salome/gui/GEOM/input/tui_test_all.doc

@@ -108,6 +108,9 @@
 \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")
 
 */

idl/GEOM_Gen.idl

@@ -4716,6 +4716,24 @@ module GEOM
      */
     double MinSurfaceCurvatureByPoint (in GEOM_Object theShape, in GEOM_Object thePoint);
 
+    /*!
+     *  \brief Get vector of curvature of surface in the given point along the given direction.
+     *  \param theShape - face.
+     *  \param thePoint - point.
+     *  \param theDirection - direction.
+     *  \note  Before the calculation of curvature, the point and the direction
+     *         are projected to the face, if the point does not lay on it or
+     *         the direction is not tangent to it initially.
+     *  \return Vector of curvature. The returned vector is codirectional with
+     *          the normal to the 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.
+     */
+    GEOM_Object SurfaceCurvatureByPointAndDirection (in GEOM_Object theShape,
+                                                     in GEOM_Object thePoint,
+                                                     in GEOM_Object theDirection);
+
   };
 
  // # GEOM_IGroupOperations:

src/GEOMGUI/GEOM_msg_en.ts

@@ -4980,6 +4980,10 @@ Please, select face, shell or solid and try again</translation>
         <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>

src/GEOMGUI/GEOM_msg_fr.ts

@@ -4972,6 +4972,10 @@ Choisissez une face, une coque ou un solide et essayez de nouveau</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>

src/GEOMGUI/GEOM_msg_ja.ts

@@ -4975,6 +4975,10 @@
       <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>

src/GEOMImpl/GEOMImpl_IMeasure.hxx

@@ -30,7 +30,8 @@ class GEOMImpl_IMeasure
     MEASURE_ARG_BASE  = 1,
     MEASURE_ARG_POINT = 2,
     MEASURE_INDEX = 3,
-    MEASURE_USE_ORI = 4
+    MEASURE_USE_ORI = 4,
+    MEASURE_ARG_DIR = 5
   };
  public:
 
@@ -56,6 +57,11 @@ class GEOMImpl_IMeasure
                               !_func->IsDone() ); // old behavior was to useOri
   }
 
+  void SetDirection(Handle(GEOM_Function) theDir)
+  { _func->SetReference(MEASURE_ARG_DIR, theDir); }
+
+  Handle(GEOM_Function) GetDirection() { return _func->GetReference(MEASURE_ARG_DIR); }
+
  private:
 
   Handle(GEOM_Function) _func;

src/GEOMImpl/GEOMImpl_IMeasureOperations.cxx

@@ -2658,6 +2658,63 @@ Standard_Real GEOMImpl_IMeasureOperations::MinSurfaceCurvatureByPoint
   return getSurfaceCurvatures(aSurf, UV.X(), UV.Y(), false);
 }
 
+//=============================================================================
+/*!
+ *  SurfaceCurvatureByPointAndDirection
+ */
+//=============================================================================
+Handle(GEOM_Object) GEOMImpl_IMeasureOperations::SurfaceCurvatureByPointAndDirection
+                                                 (Handle(GEOM_Object) theSurf,
+                                                  Handle(GEOM_Object) thePoint,
+                                                  Handle(GEOM_Object) theDirection)
+{
+  SetErrorCode(KO);
+
+  if (theSurf.IsNull() || thePoint.IsNull() || theDirection.IsNull()) return NULL;
+
+  Handle(GEOM_Function) aSurf = theSurf->GetLastFunction();
+  Handle(GEOM_Function) aPoint = thePoint->GetLastFunction();
+  Handle(GEOM_Function) aDirection = theDirection->GetLastFunction();
+  if (aSurf.IsNull() || aPoint.IsNull() || aDirection.IsNull()) return NULL;
+
+  //Add a new CurvatureVector object
+  //Handle(GEOM_Object) aCV = GetEngine()->AddObject(GEOM_CURVATURE_VEC);
+  Handle(GEOM_Object) aCV = GetEngine()->AddObject(GEOM_VECTOR);
+
+  //Add a new CurvatureVector function
+  Handle(GEOM_Function) aFunction =
+    aCV->AddFunction(GEOMImpl_MeasureDriver::GetID(), CURVATURE_VEC_MEASURE);
+  if (aFunction.IsNull()) return NULL;
+
+  //Check if the function is set correctly
+  if (aFunction->GetDriverGUID() != GEOMImpl_MeasureDriver::GetID()) return NULL;
+
+  GEOMImpl_IMeasure aCI (aFunction);
+  aCI.SetBase(aSurf);
+  aCI.SetPoint(aPoint);
+  aCI.SetDirection(aDirection);
+
+  //Compute the CurvatureVector
+  try {
+    OCC_CATCH_SIGNALS;
+    if (!GetSolver()->ComputeFunction(aFunction)) {
+      SetErrorCode("Measure driver failed to compute a surface curvature");
+      return NULL;
+    }
+  }
+  catch (Standard_Failure& aFail) {
+    SetErrorCode(aFail.GetMessageString());
+    return NULL;
+  }
+
+  //Make a Python command
+  GEOM::TPythonDump(aFunction) << aCV << " = geompy.CurvatureOnFace(" << theSurf
+                               << ", " << thePoint << ", " << theDirection << ")";
+
+  SetErrorCode(OK);
+  return aCV;
+}
+
 //=======================================================================
 //function : FillErrorsSub
 //purpose  : Fill the errors list of subshapes on shape.

src/GEOMImpl/GEOMImpl_IMeasureOperations.hxx

@@ -216,6 +216,11 @@ class GEOMImpl_IMeasureOperations : public GEOM_IOperations {
   Standard_EXPORT Standard_Real MinSurfaceCurvatureByPoint (Handle(GEOM_Object) theSurf,
                                                             Handle(GEOM_Object) thePoint);
 
+  Standard_EXPORT Handle(GEOM_Object) SurfaceCurvatureByPointAndDirection
+                                      (Handle(GEOM_Object) theSurf,
+                                       Handle(GEOM_Object) thePoint,
+                                       Handle(GEOM_Object) theDirection);
+
  private:
 
    void FillErrorsSub

src/GEOMImpl/GEOMImpl_MeasureDriver.cxx

@@ -33,10 +33,12 @@
 #include <BRep_Tool.hxx>
 #include <BRepTools.hxx>
 #include <BRepGProp.hxx>
+#include <BRepLProp_SLProps.hxx>
 #include <BRepBuilderAPI_Copy.hxx>
 #include <BRepBuilderAPI_MakeEdge.hxx>
 #include <BRepBuilderAPI_MakeVertex.hxx>
 #include <BRepPrimAPI_MakeBox.hxx>
+#include <BRepOffsetAPI_NormalProjection.hxx>
 
 #include <TopAbs.hxx>
 #include <TopoDS.hxx>
@@ -49,7 +51,9 @@
 
 #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>
@@ -80,6 +84,91 @@ GEOMImpl_MeasureDriver::GEOMImpl_MeasureDriver()
 {
 }
 
+//! This function is designed to evaluate normal curvature of the surface
+//! in the given point along the given direction.
+//! param[in] theFace face of interest.
+//! param[in] thePoint point of interest.
+//! param[in] theDir edge, giving the direction of interest.
+//! return Edge, representing the curvature vector
+TopoDS_Shape EvaluateAlongCurvature(const TopoDS_Shape& theFace,
+                                    const TopoDS_Shape& thePoint,
+                                    const TopoDS_Shape& theDir)
+{
+  // Point
+  if (thePoint.IsNull())
+    Standard_NullObject::Raise("Point for curvature measurement is null");
+  if (thePoint.ShapeType() != TopAbs_VERTEX)
+    Standard_TypeMismatch::Raise("Point for curvature calculation is not a vertex");
+  gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(thePoint));
+
+  // Point projection on the face
+  Standard_Real U, V;
+  aPnt = GEOMUtils::ProjectPointOnFace(aPnt, theFace, U, V);
+  gp_Pnt2d UV (U, V);
+
+  // 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.IsNormalDefined())
+    Standard_ConstructionError::Raise
+      ("Curvature calculation failed: normal direction is not defined");
+
+  // Get differential properties
+  gp_Vec n = Props.Normal();
+  if (aV.IsParallel(n, Precision::Angular()))
+    Standard_ConstructionError::Raise
+      ("Curvature calculation failed: direction is normal to the face");
+
+  if (!Props.IsCurvatureDefined())
+    Standard_ConstructionError::Raise
+      ("Curvature calculation failed: curvature is not defined");
+
+  // Curvature
+  Standard_Real k = 0.;
+
+  if (Props.IsUmbilic()) {
+    k = Props.MaxCurvature();
+  }
+  else {
+    gp_Dir aMaxDir, aMinDir;
+    Props.CurvatureDirections(aMaxDir, aMinDir);
+
+    gp_Vec2d T (aV.Dot(aMinDir), aV.Dot(aMaxDir));
+    if (Abs(T.X()) < Precision::Confusion() &&
+        Abs(T.Y()) < Precision::Confusion())
+      Standard_ConstructionError::Raise
+                     ("Curvature calculation failed: direction is normal to the face");
+    gp_Dir2d aDirT (T);
+    Standard_Real cosAlpha = aDirT.X();
+
+    Standard_Real aMinCurv = Props.MinCurvature();
+    Standard_Real aMaxCurv = Props.MaxCurvature();
+
+    // Euler formula: k = k1 * cos^2(alpha) + k2 * sin^2(alpha)
+    //                k = k2 + (k1 - k2) * cos^2(alpha)
+    k = aMaxCurv + (aMinCurv - aMaxCurv) * cosAlpha * cosAlpha;
+  }
+
+  if (Abs(k) < Precision::Confusion())
+    Standard_Failure::Raise("ZERO_CURVATURE");
+
+  // Radius or -radius of curvature
+  Standard_Real r = 1.0 / k;
+
+  // Result
+  gp_Dir aNormal (n);
+  gp_Pnt aPntEnd (aPnt.XYZ() + aNormal.XYZ() * r);
+  BRepBuilderAPI_MakeEdge aBuilder (aPnt, aPntEnd);
+  if (!aBuilder.IsDone())
+    Standard_ConstructionError::Raise("Curvature calculation failed: edge is not built");
+
+  return aBuilder.Shape();
+}
+
 //=======================================================================
 //function : Execute
 //purpose  :
@@ -310,6 +399,17 @@ Standard_Integer GEOMImpl_MeasureDriver::Execute(Handle(TFunction_Logbook)& log)
       Standard_NullObject::Raise("Vector construction failed");
     aShape = aBuilder.Shape();
   }
+  else if (aType == CURVATURE_VEC_MEASURE) {
+    Handle(GEOM_Function) aSrfFunc = aCI.GetBase();
+    Handle(GEOM_Function) aPntFunc = aCI.GetPoint();
+    Handle(GEOM_Function) aDirFunc = aCI.GetDirection();
+
+    TopoDS_Shape aFace   = aSrfFunc->GetValue();
+    TopoDS_Shape aVertex = aPntFunc->GetValue();
+    TopoDS_Shape anEdge  = aDirFunc->GetValue();
+
+    aShape = EvaluateAlongCurvature(aFace, aVertex, anEdge);
+  }
   else {
   }
 
@@ -359,6 +459,12 @@ GetCreationInformation(std::string&             theOperationName,
     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;
   }

src/GEOMImpl/GEOMImpl_ProjectionDriver.cxx

@@ -127,70 +127,11 @@ Standard_Integer GEOMImpl_ProjectionDriver::Execute(Handle(TFunction_Logbook)& l
     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_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(anOriginal));
-      gp_Pnt2d aProjPnt (U, V);
+      gp_Pnt surfPnt = GEOMUtils::ProjectPointOnFace(aPnt, aFaceShape, 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);
-
       aShape = BRepBuilderAPI_MakeVertex(surfPnt).Shape();
     }
     else {

src/GEOMImpl/GEOMImpl_Types.hxx

@@ -121,6 +121,8 @@
 
 #define GEOM_EXTRACTION 58
 
+#define GEOM_CURVATURE_VEC 59
+
 //GEOM_Function types
 
 #define COPY_WITH_REF    1
@@ -357,6 +359,7 @@
 #define BND_BOX_MEASURE_PRECISE 3
 #define VECTOR_FACE_NORMALE 4
 #define VERTEX_BY_INDEX 5
+#define CURVATURE_VEC_MEASURE 6
 
 #define GROUP_FUNCTION 1
 

src/GEOMUtils/GEOMUtils.cxx

@@ -37,6 +37,7 @@
 #include <BRepGProp.hxx>
 #include <BRepTools.hxx>
 
+#include <BRepClass_FaceClassifier.hxx>
 #include <BRepClass3d_SolidClassifier.hxx>
 
 #include <BRepBuilderAPI_MakeFace.hxx>
@@ -64,6 +65,8 @@
 #include <TopTools_ListIteratorOfListOfShape.hxx>
 #include <TopTools_Array1OfShape.hxx>
 
+#include <GeomAPI_ProjectPointOnSurf.hxx>
+
 #include <Geom_Circle.hxx>
 #include <Geom_Surface.hxx>
 #include <Geom_Plane.hxx>
@@ -1021,6 +1024,65 @@ Standard_Real GEOMUtils::GetMinDistance
   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

src/GEOMUtils/GEOMUtils.hxx

@@ -215,6 +215,19 @@ namespace GEOMUtils
 					       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.
    *

src/GEOM_I/GEOM_IMeasureOperations_i.cc

@@ -1188,3 +1188,32 @@ CORBA::Double GEOM_IMeasureOperations_i::MinSurfaceCurvatureByPoint
 
   return GetOperations()->MinSurfaceCurvatureByPoint(aShape,aPoint);
 }
+
+//=============================================================================
+/*!
+ *  SurfaceCurvatureByPointAndDirection
+ */
+//=============================================================================
+GEOM::GEOM_Object_ptr GEOM_IMeasureOperations_i::SurfaceCurvatureByPointAndDirection
+                                          (GEOM::GEOM_Object_ptr theSurf,
+                                           GEOM::GEOM_Object_ptr thePoint,
+                                           GEOM::GEOM_Object_ptr theDirection)
+{
+  GEOM::GEOM_Object_var aGEOMObject;
+
+  //Set a not done flag
+  GetOperations()->SetNotDone();
+
+  //Get the reference shape
+  Handle(::GEOM_Object) aShape = GetObjectImpl(theSurf);
+  Handle(::GEOM_Object) aPoint = GetObjectImpl(thePoint);
+  Handle(::GEOM_Object) aDirection = GetObjectImpl(theDirection);
+  if (aShape.IsNull() || aPoint.IsNull() || aDirection.IsNull()) return aGEOMObject._retn();
+
+  Handle(::GEOM_Object) anObject =
+    GetOperations()->SurfaceCurvatureByPointAndDirection(aShape,aPoint,aDirection);
+  if (!GetOperations()->IsDone() || anObject.IsNull())
+    return aGEOMObject._retn();
+
+  return GetObject(anObject);
+}

src/GEOM_I/GEOM_IMeasureOperations_i.hh

@@ -164,6 +164,10 @@ class GEOM_I_EXPORT GEOM_IMeasureOperations_i :
   CORBA::Double MinSurfaceCurvatureByPoint (GEOM::GEOM_Object_ptr theSurf,
                                             GEOM::GEOM_Object_ptr thePoint);
 
+  GEOM::GEOM_Object_ptr SurfaceCurvatureByPointAndDirection (GEOM::GEOM_Object_ptr theSurf,
+                                                             GEOM::GEOM_Object_ptr thePoint,
+                                                             GEOM::GEOM_Object_ptr theDirection);
+
   ::GEOMImpl_IMeasureOperations* GetOperations()
   { return (::GEOMImpl_IMeasureOperations*)GetImpl(); }
 };

src/GEOM_SWIG/GEOM_TestAll.py

@@ -595,5 +595,20 @@ def TestAll (geompy, math):
   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)
+  curvature_zero = geompy.CurvatureOnFace(Face_2, px, vz)
+  assert(geompy.MeasuOp.GetErrorCode() == "ZERO_CURVATURE")
+  assert(not curvature_zero)
+  isExcept = False
+  try:
+    # p0 is on cylinder axis, projection should fail
+    geompy.CurvatureOnFace(Face_2, p0, vy)
+  except:
+    isExcept = True
+  assert(isExcept)
 
   print("DONE")

src/GEOM_SWIG/geomBuilder.py

@@ -11178,6 +11178,56 @@ class geomBuilder(GEOM._objref_GEOM_Gen):
             return aSurf
         ## @}
 
+        ## 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.
+        #  @param theName Object name; when specified, this parameter is used
+        #         for result publication in the study. Otherwise, if automatic
+        #         publication is switched on, default value is used for result name.
+        #
+        #  @return New GEOM.GEOM_Object, containing vector of curvature of theSurf.
+        #          The returned vector is codirectional with the normal to the 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 radius.
+        #          Null shape is returned in case of infinite radius
+        #          (zero curvature), for example, in case of flat face.
+        #
+        ## @ref swig_CurvatureOnFace "Example"
+        @ManageTransactions("MeasuOp")
+        def CurvatureOnFace(self, theSurf, thePoint, theDirection, theName=None):
+            """
+            Measure curvature radius of surface in the given point along the given direction.
+
+            Parameters:
+                theSurf the given face.
+                thePoint given point.
+                theDirection given direction.
+                theName Object name; when specified, this parameter is used
+                        for result publication in the study. Otherwise, if automatic
+                        publication is switched on, default value is used for result name.
+
+            Returns:
+                New GEOM.GEOM_Object, containing vector of curvature of theSurf.
+                The returned vector is codirectional with the normal to the 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 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)
+            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
         #  @param theShape Shape, to get tolerances of.
         #  @return [FaceMin,FaceMax, EdgeMin,EdgeMax, VertMin,VertMax]\n