[SofaImplicitField] Add getHessian and getGradient to the python binding (#5655)

* Add binding for hessian & gradient.

* Add getHessian, getGradient and getValue with an unified API to Binding_ScalarField

I unified the interface so it is conformant with the one from sofa expect for getHessian.

* Update applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp

Author: damienmarchal

applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py

@@ -1,5 +1,6 @@
 import Sofa
 from SofaImplicitField import ScalarField
+from SofaTypes.SofaTypes import Vec3d, Mat3x3
 import numpy
 
 class Sphere(ScalarField):
@@ -9,20 +10,41 @@ def __init__(self, *args, **kwargs):
         self.addData("center", type="Vec3d",value=kwargs.get("center", [0.0,0.0,0.0]), default=[0.0,0.0,0.0], help="center of the sphere", group="Geometry")
         self.addData("radius", type="double",value=kwargs.get("radius", 1.0), default=1, help="radius of the sphere", group="Geometry")
 
-    def getValue(self, x, y, z):
+    def getValue(self, position):
+        x,y,z = position
         return numpy.sqrt( numpy.sum((self.center.value - numpy.array([x,y,z]))**2) ) - self.radius.value 
 
+class SphereWithCustomHessianAndGradient(Sphere):
+    def __init__(self, *args, **kwargs):
+        Sphere.__init__(self, *args, **kwargs)#
+
+    def getGradient(self, position):
+        return Vec3d(3.0,2.0,1.0)
+
+    def getHessian(self, position):
+        return Mat3x3([[1,1,1],[2,1,1],[3,1,1]])
+
 class FieldController(Sofa.Core.Controller):
     def __init__(self, *args, **kwargs):
         Sofa.Core.Controller.__init__(self, *args, **kwargs)
         self.field = kwargs.get("target")
 
     def onAnimateEndEvent(self, event):
-        print("Animation end event")
-        print("Field value at 0,0,0 is: ", self.field.getValue(0.0,0.0,0.0) ) 
-        print("Field value at 1,0,0 is: ", self.field.getValue(1.0,0.0,0.0) ) 
-        print("Field value at 2,0,0 is: ", self.field.getValue(2.0,0.0,0.0) ) 
+        print("Animation end event, ")
+        print("Field value at 0,0,0 is: ", self.field.getValue(Vec3d(0.0,0.0,0.0)) )
+        print("Field value at 1,0,0 is: ", self.field.getValue(Vec3d(1.0,0.0,0.0)) )
+        print("Field value at 2,0,0 is: ", self.field.getValue(Vec3d(2.0,0.0,0.0)) )
+
+        print("Gradient value at 0,0,0 is: ", type(self.field.getGradient(Vec3d(0.0,0.0,0.0))))
+        print("Hessian value at 0,0,0 is: ", type(self.field.getHessian(Vec3d(0.0,0.0,0.0))))
 
 def createScene(root):
-    root.addObject(Sphere("field"))  
-    root.addObject(FieldController(target=root.field))
+    """In this scene we create two scalar field of spherical shape, the two are implemented using
+    python. The first one is overriding only the getValue, the hessian and gradient is thus computed using
+    finite difference in the c++ code. The second field is overriding the hessian and gradient function
+    """
+    root.addObject(Sphere("field1"))
+    root.addObject(FieldController(name="controller1", target=root.field1))
+
+    root.addObject(SphereWithCustomHessianAndGradient("field2"))
+    root.addObject(FieldController(name="controller2", target=root.field2))

applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp

@@ -33,22 +33,65 @@ namespace sofaimplicitfield {
 using namespace sofapython3;
 using sofa::component::geometry::ScalarField;
 using sofa::core::objectmodel::BaseObject;
-using sofa::type::Vec3d;
+using sofa::type::Vec3;
+using sofa::type::Mat3x3;
 
 class ScalarField_Trampoline : public ScalarField {
 public:
     SOFA_CLASS(ScalarField_Trampoline, ScalarField);
 
-    double getValue(Vec3d& pos, int& domain) override{
+    // Override this function so that it returns the actual python class name instead of
+    // "ScalarField_Trampoline" which correspond to this utility class.
+    std::string getClassName() const override
+    {
+        PythonEnvironment::gil acquire;
+
+        // Get the actual class name from python.
+        return py::str(py::cast(this).get_type().attr("__name__"));
+    }
+
+    double getValue(Vec3& pos, int& domain) override
+    {
+        SOFA_UNUSED(domain);
+        PythonEnvironment::gil acquire;
+
+        PYBIND11_OVERLOAD_PURE(double, ScalarField, getValue, pos);
+    }
+
+    Vec3 getGradient(Vec3& pos, int& domain) override
+    {
         SOFA_UNUSED(domain);
         PythonEnvironment::gil acquire;
 
-        PYBIND11_OVERLOAD_PURE(double, ScalarField, getValue, pos.x(), pos.y(), pos.z());
+        PYBIND11_OVERLOAD(Vec3, ScalarField, getGradient, pos);
+    }
+
+    void getHessian(Vec3 &pos, Mat3x3& h) override
+    {
+        /// The implementation is a bit more complex compared to getGradient. This is because we change de signature between the c++ API and the python one.
+        PythonEnvironment::gil acquire;
+
+        // Search if there is a python override,
+        pybind11::function override = pybind11::get_override(static_cast<const ScalarField*>(this),"getHessian");
+        if(!override){
+            return ScalarField::getHessian(pos, h);
+        }
+        // as there is one override, we call it, passing the "pos" argument and storing the return of the
+        // value in the "o" variable.
+        auto o = override(pos);
+
+        // then we check that the function correctly returned a Mat3x3 object and copy its value
+        // in case there is no Mat3x3 returned values, rise an error
+        if(py::isinstance<Mat3x3>(o))
+            h = py::cast<Mat3x3>(o);
+        else
+            throw py::type_error("The function getHessian must return a Mat3x3");
+        return;
     }
 };
 
 void moduleAddScalarField(py::module &m) {
-    py::class_<ScalarField, BaseObject, ScalarField_Trampoline,
+    py::class_<ScalarField, ScalarField_Trampoline, BaseObject,
                py_shared_ptr<ScalarField>> f(m, "ScalarField", py::dynamic_attr(), "");
 
     f.def(py::init([](py::args &args, py::kwargs &kwargs) {
@@ -74,7 +117,23 @@ void moduleAddScalarField(py::module &m) {
         return ff;
     }));
 
-    m.def("getValue", &ScalarField_Trampoline::getValue);
+    f.def("getValue", [](ScalarField* self, Vec3 pos){
+        int domain=-1;
+        // This shouldn't be self->ScalarField::getValue because it is a pure function
+        // so there is not ScalarField::getValue emitted.
+        return self->getValue(pos, domain);
+    });
+
+    f.def("getGradient", [](ScalarField* self, Vec3 pos){
+        int domain=-1;
+        return self->ScalarField::getGradient(pos, domain);
+    });
+
+    f.def("getHessian", [](ScalarField* self, Vec3 pos){
+        Mat3x3 result;
+        self->getHessian(pos, result);
+        return result;
+    });
 }
 
 }