Add vectorized version of the getValue function in ScalarField.

examples/python/xshape/operators.py

@@ -12,6 +12,10 @@ def __init__(self, *args, **kwargs):
     def getValue(self, position):
         return min(self.childA.getValue(position), self.childB.getValue(position))
 
+    def getValues(self, positions, results):
+        results[:] = numpy.minimum(self.childA.getValues(positions, numpy.empty(len(positions))), self.childB.getValues(positions, numpy.empty(len(positions))))
+        return results
+
 class Difference(ScalarField):
     """Difference of two scalar fields"""
     def __init__(self, *args, **kwargs):
@@ -23,6 +27,10 @@ def __init__(self, *args, **kwargs):
     def getValue(self, position):
         return max(-self.childA.getValue(position), self.childB.getValue(position))
 
+    def getValues(self, positions, results):
+        results[:] = numpy.maximum(-self.childA.getValues(positions, numpy.empty(len(positions))), self.childB.getValues(positions, numpy.empty(len(positions))))
+        return results
+
 class Intersection(ScalarField):
     """Intersection of two scalar fields"""
     def __init__(self, *args, **kwargs):
@@ -32,4 +40,8 @@ def __init__(self, *args, **kwargs):
         self.childB = kwargs.get("childB", None)
 
     def getValue(self, position):
-        return max(self.childA.getValue(position), self.childB.getValue(position))
+        return max(self.childA.getValue(position), self.childB.getValue(position))
+
+    def getValues(self, positions, results):
+        results[:] = numpy.maximum(self.childA.getValues(positions, numpy.empty(len(positions))), self.childB.getValues(positions, numpy.empty(len(positions))))
+        return results  

examples/python/xshape/primitives.py

@@ -14,9 +14,23 @@ def __init__(self, *args, **kwargs):
         self.addData("radius", type="double",value=kwargs.get("radius", 1.0), default=1, help="radius of the sphere", group="Geometry")
 
     def getValue(self, pos):
+        """This version is of very low performance as there are a huge amount of call to the python side"""
         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
+
+    def getValues(self, positions, results):
+        """This version of the overrides the getValues so that we fetch the data once"""
+        results[:] = numpy.linalg.norm(positions - self.center.value, axis=1) - self.radius.value
+        return results
+
 class RoundedBox(ScalarField):
     def __init__(self, *args, **kwargs):
         ScalarField.__init__(self, *args, **kwargs)
@@ -32,3 +46,12 @@ def getValue(self, pos):
         q = numpy.abs(self.center.value - numpy.array([x,y,z])) - b + r
         res = numpy.linalg.norm(numpy.maximum(q, 0.0)) + min(max(q[0], max(q[1],q[2]) ), 0.0) - r
         return res
+
+    def getValues(self, positions, results):
+        b = self.dimensions.value
+        r = self.rounding_radius.value
+        q = numpy.abs(positions - self.center.value) - b + r
+        outside = numpy.linalg.norm(numpy.maximum(q, 0.0), axis=1)
+        inside = numpy.minimum(numpy.max(q, axis=1), 0.0)
+        results[:] = outside + inside - r
+        return results

python/src/Binding_ScalarField.cpp

@@ -36,6 +36,59 @@ 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 +111,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);
@@ -125,6 +196,10 @@ void moduleAddScalarField(py::module &m) {
         return self->getValue(pos, domain);
     });
 
+    f.def("getValues", [](ScalarField* self, const std::vector<Vec3>& positions, std::vector<double>& results){
+        self->ScalarField::getValues(positions, results);
+    });
+
     f.def("getGradient", [](ScalarField* self, Vec3 pos){
         int domain=-1;
         return self->ScalarField::getGradient(pos, domain);

src/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;
 }
 
@@ -110,11 +111,7 @@ void FieldToSurfaceMesh::updateMeshIfNeeded()
 
     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->getValues(positions, res);
                                       },
                                      tmpPoints, tmpTriangles);
 

src/SofaImplicitField/components/engine/FieldToSurfaceMesh.h

@@ -35,15 +35,16 @@ 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 ;

src/SofaImplicitField/components/geometry/ScalarField.cpp

@@ -68,6 +68,16 @@ Vec3d ScalarField::getGradientByFinitDifference(Vec3d& pos, int& i)
     return Result;
 }
 
+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::getGradient(Vec3d& pos, int& i)
 {
     return getGradientByFinitDifference(pos, i);

src/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);