Add vectorized version

Author: damienmarchal

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/geometry/ScalarField.cpp

@@ -45,6 +45,17 @@ 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));
+    }
+    return;
+}
+
 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/xshape/primitives.py

@@ -5,6 +5,7 @@
 """
 from SofaImplicitField import ScalarField
 import numpy
+import numba
 
 class Sphere(ScalarField):
     def __init__(self, *args, **kwargs):
@@ -17,6 +18,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);