diff --git a/examples/python/xshape/operators.py b/examples/python/xshape/operators.py index 4506c13..ea7dda5 100644 --- a/examples/python/xshape/operators.py +++ b/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)) \ No newline at end of file + 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 \ No newline at end of file diff --git a/examples/python/xshape/primitives.py b/examples/python/xshape/primitives.py index 87f8a34..b4d6422 100644 --- a/examples/python/xshape/primitives.py +++ b/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 \ No newline at end of file diff --git a/python/src/Binding_ScalarField.cpp b/python/src/Binding_ScalarField.cpp index 927a801..e9ae6d3 100644 --- a/python/src/Binding_ScalarField.cpp +++ b/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 vector_to_numpy(const std::vector& 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 shape = { size, 3 }; + + // Strides en bytes (ici, contiguous : 3 colonnes, chaque double 8 octets) + std::vector strides = { 3 * sizeof(double), sizeof(double) }; + + return py::array_t( + shape, + strides, + data_ptr, // data pointer + free_when_done // capsule pour gérer la vie mémoire + ); +} + +py::array_t vector_to_numpy(const std::vector& 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 shape = { size }; + + // Strides en bytes (ici, contiguous : 3 colonnes, chaque double 8 octets) + std::vector strides = { sizeof(double) }; + + return py::array_t( + 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& positions, std::vector& results) override + { + PythonEnvironment::gil acquire; + + // Search if there is a python override, + pybind11::function override = pybind11::get_override(static_cast(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& positions, std::vector& results){ + self->ScalarField::getValues(positions, results); + }); + f.def("getGradient", [](ScalarField* self, Vec3 pos){ int domain=-1; return self->ScalarField::getGradient(pos, domain); diff --git a/src/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp b/src/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp index e42e184..4864cd6 100644 --- a/src/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp +++ b/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& positions, std::vector& res){ - int i=0; - for(auto& position : positions) - { - res[i++]=field->getValue(position); - } + field->getValues(positions, res); }, tmpPoints, tmpTriangles); diff --git a/src/SofaImplicitField/components/engine/FieldToSurfaceMesh.h b/src/SofaImplicitField/components/engine/FieldToSurfaceMesh.h index e92783b..69c03ba 100644 --- a/src/SofaImplicitField/components/engine/FieldToSurfaceMesh.h +++ b/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 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 ; diff --git a/src/SofaImplicitField/components/geometry/ScalarField.cpp b/src/SofaImplicitField/components/geometry/ScalarField.cpp index 9dcae38..de0396a 100644 --- a/src/SofaImplicitField/components/geometry/ScalarField.cpp +++ b/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& positions, std::vector& 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); diff --git a/src/SofaImplicitField/components/geometry/ScalarField.h b/src/SofaImplicitField/components/geometry/ScalarField.h index d676dd2..d603d2d 100644 --- a/src/SofaImplicitField/components/geometry/ScalarField.h +++ b/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& positions, std::vector& 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);