[SofaImplicitField] add vectorized scalarfield binding

applications/plugins/SofaImplicitField/MarchingCube.cpp

@@ -43,7 +43,6 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep,
     if( nz < 2 || ny < 2 || nx < 2 )
         return;
 
-    double cx,cy,cz;
     int z,mk;
     const int *tri;
 

applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp

@@ -66,6 +66,7 @@ void FieldToSurfaceMesh::init()
         d_componentState = core::objectmodel::ComponentState::Invalid;
     }
 
+    updateMeshIfNeeded();
     d_componentState = core::objectmodel::ComponentState::Valid;
 }
 
@@ -109,12 +110,8 @@ void FieldToSurfaceMesh::updateMeshIfNeeded()
     tmpTriangles.clear();
 
     marchingCube.generateSurfaceMesh(isoval, mstep, invStep, gridmin, gridmax,
-                                     [field](std::vector<Vec3d>& positions, std::vector<double>& res){
-                                        int i=0;
-                                        for(auto& position : positions)
-                                        {
-                                            res[i++]=field->getValue(position);
-                                        }
+                                     [field](const std::vector<Vec3d>& positions, std::vector<double>& res){
+                                        field->getValues(positions, res);
                                       },
                                      tmpPoints, tmpTriangles);
 

applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h

@@ -24,7 +24,6 @@
 #include <SofaImplicitField/components/geometry/ScalarField.h>
 #include <SofaImplicitField/MarchingCube.h>
 #include <sofa/core/topology/BaseMeshTopology.h>
-#include <future>
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 namespace sofaimplicitfield::component::engine
@@ -35,15 +34,15 @@ typedef sofa::core::topology::BaseMeshTopology::SeqTriangles SeqTriangles;
 typedef sofa::core::topology::BaseMeshTopology::Triangle Triangle;
 typedef sofa::type::vector<sofa::type::Vec3d> VecCoord;
 
+using sofa::core::objectmodel::BaseComponent;
 using sofa::component::geometry::ScalarField;
 using sofa::core::visual::VisualParams ;
-using BaseObject [[deprecated("Use sofa::core::objectmodel::BaseObject instead.")]] = sofa::core::objectmodel::BaseObject;
 using sofa::type::Vec3d ;
 
-class FieldToSurfaceMesh : public BaseObject
+class FieldToSurfaceMesh : public BaseComponent
 {
 public:
-    SOFA_CLASS(FieldToSurfaceMesh, BaseObject);
+    SOFA_CLASS(FieldToSurfaceMesh, BaseComponent);
 
     virtual void init() override ;
     virtual void draw(const VisualParams*params) override ;

applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp

@@ -45,6 +45,16 @@ void ScalarField::init()
     d_componentState.setValue(core::objectmodel::ComponentState::Valid);
 }
 
+void ScalarField::getValues(const std::vector<Vec3d>& positions, std::vector<double>& results)
+{
+    results.clear();
+    results.reserve(positions.size());
+    for(auto position : positions)
+    {
+        results.emplace_back(getValue(position));
+    }
+}
+
 Vec3d ScalarField::getGradientByFinitDifference(Vec3d& pos, int& i)
 {
     Vec3d Result;

applications/plugins/SofaImplicitField/components/geometry/ScalarField.h

@@ -66,6 +66,9 @@ class SOFA_SOFAIMPLICITFIELD_API ScalarField : public BaseObject
     virtual double getValue(Vec3d& pos, int& domain) = 0;
     inline double getValue(Vec3d& pos) { int domain=-1; return getValue(pos,domain); }
 
+    // Compute the field for a range or input values
+    virtual void getValues(const std::vector<Vec3d>& positions, std::vector<double>& results);
+
     /// By default compute the gradient using a first order finite difference approache
     /// If you have analytical derivative don't hesitate to override this function.
     virtual Vec3d getGradient(Vec3d& pos, int& domain);

applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py

@@ -18,8 +18,7 @@ def createScene(root : Sofa.Core.Node):
        The scalar fields are 'spherical', one implemented in python, the other in c++
        One of the produced mesh is then connected to a visual model.
     """
-    root.addObject("RequiredPlugin", pluginName="SofaImplicitField")
-
+    root.addObject("RequiredPlugin", name="SofaImplicitField")
     root.addObject(DrawController())
 
     ########################### Fields ##################

applications/plugins/SofaImplicitField/examples/python/python-implicit-field-example.py

@@ -0,0 +1,17 @@
+import Sofa
+from primitives import Sphere
+
+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]) )
+
+def createScene(root):
+    root.addObject(Sphere("field"))
+    root.addObject(FieldController(target=root.field))

applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py

@@ -17,6 +17,14 @@ def getValue(self, pos):
         x,y,z = pos
         return numpy.linalg.norm(self.center.value - numpy.array([x,y,z])) - self.radius.value
 
+    def getValues(self, positions, out_values):
+        """This version of the overrides the getValues so that we fetch the data once"""        
+        center = self.center.value
+        radius = self.radius.value
+        for i in range(len(positions)):
+            r = numpy.linalg.norm(center - positions[i]) - radius
+            out_values[i] = r
+
 class RoundedBox(ScalarField):
     def __init__(self, *args, **kwargs):
         ScalarField.__init__(self, *args, **kwargs)

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

@@ -36,6 +36,58 @@ using sofa::core::objectmodel::BaseObject;
 using sofa::type::Vec3;
 using sofa::type::Mat3x3;
 
+py::array_t<double> vector_to_numpy(const std::vector<Vec3>& vec) {
+    // Pybind11 gère la mémoire en créant un capsule pour que numpy sache comment libérer
+    // On transfère la propriété du vecteur à numpy, donc pas de copie
+    const double* data_ptr = (const double*)vec.data();
+    size_t size = vec.size();
+
+    // capsule: mémoire gérée par le vector, sera libérée quand python détruit l'objet numpy
+    py::capsule free_when_done(vec.data(), [](void *) {
+        // On ne fait rien ici car vector gère sa mémoire.
+        // Si on voulait transférer la propriété, on ferait delete ici.
+    });
+
+    // Dimensions du tableau numpy
+    std::vector<size_t> shape = { size, 3 };
+
+    // Strides en bytes (ici, contiguous : 3 colonnes, chaque double 8 octets)
+    std::vector<size_t> strides = { 3 * sizeof(double), sizeof(double) };
+
+    return py::array_t<double>(
+        shape,
+        strides,
+        data_ptr,      // data pointer
+        free_when_done // capsule pour gérer la vie mémoire
+    );
+}
+
+py::array_t<double> vector_to_numpy(const std::vector<double>& vec) {
+    // Pybind11 gère la mémoire en créant un capsule pour que numpy sache comment libérer
+    // On transfère la propriété du vecteur à numpy, donc pas de copie
+    const double* data_ptr = (const double*)vec.data();
+    size_t size = vec.size();
+
+    // capsule: mémoire gérée par le vector, sera libérée quand python détruit l'objet numpy
+    py::capsule free_when_done(vec.data(), [](void *) {
+        // On ne fait rien ici car vector gère sa mémoire.
+        // Si on voulait transférer la propriété, on ferait delete ici.
+    });
+
+    // Dimensions du tableau numpy
+    std::vector<size_t> shape = { size };
+
+    // Strides en bytes (ici, contiguous : 3 colonnes, chaque double 8 octets)
+    std::vector<size_t> strides = { sizeof(double) };
+
+    return py::array_t<double>(
+        shape,
+        strides,
+        data_ptr,      // data pointer
+        free_when_done // capsule pour gérer la vie mémoire
+    );
+}
+
 class ScalarField_Trampoline : public ScalarField {
 public:
     SOFA_CLASS(ScalarField_Trampoline, ScalarField);
@@ -58,6 +110,24 @@ class ScalarField_Trampoline : public ScalarField {
         PYBIND11_OVERLOAD_PURE(double, ScalarField, getValue, pos);
     }
 
+    void getValues(const std::vector<Vec3>& positions, std::vector<double>& results) override
+    {
+        PythonEnvironment::gil acquire;
+
+        // Search if there is a python override,
+        pybind11::function override = pybind11::get_override(static_cast<const ScalarField*>(this),"getValues");
+        if(!override){
+            return ScalarField::getValues(positions, results);
+        }
+
+        // Be sure there is enough space to hold the results
+        results.resize(positions.size());
+
+        // 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(vector_to_numpy(positions), vector_to_numpy(results));
+    }
+
     Vec3 getGradient(Vec3& pos, int& domain) override
     {
         SOFA_UNUSED(domain);