diff --git a/applications/plugins/SofaImplicitField/CMakeLists.txt b/applications/plugins/SofaImplicitField/CMakeLists.txt index 35d08e7f7b1..aca4a06b47b 100644 --- a/applications/plugins/SofaImplicitField/CMakeLists.txt +++ b/applications/plugins/SofaImplicitField/CMakeLists.txt @@ -6,12 +6,14 @@ sofa_find_package(Sofa.Component.Topology.Container.Constant REQUIRED) set(HEADER_FILES config.h.in initSofaImplicitField.h + MarchingCube.h # This is backward compatibility deprecated/SphereSurface.h deprecated/ImplicitSurfaceContainer.h # This is a backward compatibility file toward ScalarField deprecated/InterpolatedImplicitSurface.h # This is a backward compatibility file toward DiscreteGridField + components/engine/FieldToSurfaceMesh.h components/geometry/BottleField.h components/geometry/DiscreteGridField.h components/geometry/SphericalField.h @@ -23,11 +25,13 @@ set(HEADER_FILES set(SOURCE_FILES initSofaImplicitField.cpp + MarchingCube.cpp ## This is a backward compatibility.. deprecated/SphereSurface.cpp deprecated/InterpolatedImplicitSurface.cpp + components/engine/FieldToSurfaceMesh.cpp components/geometry/BottleField.cpp components/geometry/ScalarField.cpp components/geometry/DiscreteGridField.cpp diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp new file mode 100644 index 00000000000..a9b74985655 --- /dev/null +++ b/applications/plugins/SofaImplicitField/MarchingCube.cpp @@ -0,0 +1,164 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture, development version * +* (c) 2006-2025 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include + +#include +#include +#include +#include +#include + +namespace sofaimplicitfield +{ + +void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, const double invStep, + const Vec3d& gridmin, const Vec3d& gridmax, + std::function&, std::vector&)> getFieldValueAt, + SeqCoord& tmpPoints, SeqTriangles& tmpTriangles) +{ + int nx = floor((gridmax.x() - gridmin.x()) * invStep) + 1 ; + int ny = floor((gridmax.y() - gridmin.y()) * invStep) + 1 ; + int nz = floor((gridmax.z() - gridmin.z()) * invStep) + 1 ; + + // Marching cubes only works for a grid size larger than two + if( nz < 2 || ny < 2 || nx < 2 ) + return; + + int z,mk; + const int *tri; + + // Creates two planes + CubeData c{{-1,-1,-1},0}; + planes.resize(2*nx*ny); + for(size_t i=0;i &positions, std::vector& output, + double mstep, double gridmin_y, double gridmin_x, int ny, int nx, float cz, + std::vector::iterator itDestPlane) + { + for (int i=0, y = 0 ; y < ny ; ++y) + { + double cy = gridmin_y + mstep * y ; + for (int x = 0 ; x < nx ; ++x) + { + double cx = gridmin_x + mstep * x ; + positions[i++].set(cx, cy, cz ); + } + } + getFieldValueAt(positions, output) ; + + for(auto res : output){ + itDestPlane->data = res; + itDestPlane++; + } + }; + + std::vector positions; + std::vector output; + positions.resize(nx*ny); + output.resize(nx*ny); + + fillPlane(positions, output, mstep, gridmin.y(), gridmin.x(), ny, nx, gridmin.z(), P0); + for (z=1; z<=nz; ++z) + { + fillPlane(positions, output, mstep, gridmin.y(), gridmin.x(), ny, nx, gridmin.z() + mstep * z, P1); + + int edgecube[12]; + const int edgepts[12] = {0,1,0,1,0,1,0,1,2,2,2,2}; + typename std::vector::iterator base = planes.begin(); + int ip0 = P0-base; + int ip1 = P1-base; + edgecube[0] = (ip0 -dy); + edgecube[1] = (ip0 ); + edgecube[2] = (ip0 ); + edgecube[3] = (ip0-dx ); + edgecube[4] = (ip1 -dy); + edgecube[5] = (ip1 ); + edgecube[6] = (ip1 ); + edgecube[7] = (ip1-dx ); + edgecube[8] = (ip1-dx-dy); + edgecube[9] = (ip1-dy ); + edgecube[10] = (ip1 ); + edgecube[11] = (ip1-dx ); + + unsigned int di = nx; + for(int y=1; ydata>isoval)^((P1+di-dx)->data>isoval)) + { + (P1+di)->p[0] = addPoint(tmpPoints, 0, pos,gridmin, (P1+di)->data,(P1+di-dx)->data, mstep, isoval); + } + if (((P1+di)->data>isoval)^((P1+di-dy)->data>isoval)) + { + (P1+di)->p[1] = addPoint(tmpPoints, 1, pos,gridmin,(P1+di)->data,(P1+di-dy)->data, mstep, isoval); + } + if (((P1+di)->data>isoval)^((P0+di)->data>isoval)) + { + (P1+di)->p[2] = addPoint(tmpPoints, 2, pos,gridmin,(P1+di)->data,(P0+di)->data, mstep, isoval); + } + + // All points should now be created + if ((P0+di-dx-dy)->data > isoval) mk = 1; + else mk=0; + if ((P0+di -dy)->data > isoval) mk|= 2; + if ((P0+di )->data > isoval) mk|= 4; + if ((P0+di-dx )->data > isoval) mk|= 8; + if ((P1+di-dx-dy)->data > isoval) mk|= 16; + if ((P1+di -dy)->data > isoval) mk|= 32; + if ((P1+di )->data > isoval) mk|= 64; + if ((P1+di-dx )->data > isoval) mk|= 128; + + tri=sofa::helper::MarchingCubeTriTable[mk]; + while (*tri>=0) + { + typename std::vector::iterator b = base+di; + addFace(tmpTriangles, + (b+edgecube[tri[0]])->p[edgepts[tri[0]]], + (b+edgecube[tri[1]])->p[edgepts[tri[1]]], + (b+edgecube[tri[2]])->p[edgepts[tri[2]]], tmpPoints.size()); + tri+=3; + } + ++di; + } + } + std::swap(P0, P1); + } +} + +} diff --git a/applications/plugins/SofaImplicitField/MarchingCube.h b/applications/plugins/SofaImplicitField/MarchingCube.h new file mode 100644 index 00000000000..394a6163b1c --- /dev/null +++ b/applications/plugins/SofaImplicitField/MarchingCube.h @@ -0,0 +1,83 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture, development version * +* (c) 2006-2025 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +#include +#include +#include + +//////////////////////////////////////////////////////////////////////////////////////////////////// +namespace sofaimplicitfield +{ + +typedef sofa::core::topology::BaseMeshTopology::SeqTriangles SeqTriangles; +typedef sofa::core::topology::BaseMeshTopology::Triangle Triangle; +typedef sofa::type::vector SeqCoord; +using sofa::type::Vec3d; + +class MarchingCube +{ +public: + void generateSurfaceMesh(const double isoval, const double mstep, const double invStep, + const Vec3d& gridmin, const Vec3d& gridmax, + std::function &, std::vector &)> field, + SeqCoord& tmpPoints, SeqTriangles& tmpTriangles); + +private: + /// For each cube, store the vertex indices on each 3 first edges, and the data value + struct CubeData + { + int p[3]; + double data; + }; + + sofa::type::vector planes; + typename sofa::type::vector::iterator P0; /// Pointer to first plane + typename sofa::type::vector::iterator P1; /// Pointer to second plane + + int addPoint(SeqCoord& v, int i, Vec3d pos, const Vec3d& gridmin, double v0, double v1, double step, double iso) + { + pos[i] -= (iso-v0)/(v1-v0); + v.push_back( (pos * step)+gridmin ) ; + return v.size()-1; + } + + int addFace(SeqTriangles& triangles, int p1, int p2, int p3, int nbp) + { + if ((unsigned)p1<(unsigned)nbp && + (unsigned)p2<(unsigned)nbp && + (unsigned)p3<(unsigned)nbp) + { + triangles.push_back(Triangle(p1, p3, p2)); + return triangles.size()-1; + } + else + { + return -1; + } + } +}; + + +} + diff --git a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp new file mode 100644 index 00000000000..fbf12bdbf84 --- /dev/null +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp @@ -0,0 +1,171 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture, development version * +* (c) 2006-2025 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include + +#include +using sofa::core::visual::VisualParams ; + +#include +using sofa::core::RegisterObject ; + +#include + +#include "FieldToSurfaceMesh.h" + +namespace sofaimplicitfield::component::engine +{ + +FieldToSurfaceMesh::FieldToSurfaceMesh() + : l_field(initLink("field", "The scalar field to generate a mesh from.")) + , d_step(initData(&d_step,0.1,"step","Step")) + , d_IsoValue(initData(&d_IsoValue,0.0,"isoValue","Iso Value")) + , d_gridMin(initData(&d_gridMin, Vec3d(-1,-1,-1),"min","Grid Min")) + , d_gridMax(initData(&d_gridMax, Vec3d(1,1,1),"max","Grid Max")) + , d_outPoints(initData(&d_outPoints, "points", "position of the tiangles vertex")) + , d_outTriangles(initData(&d_outTriangles, "triangles", "list of triangles")) + , d_debugDraw(initData(&d_debugDraw,false, "debugDraw","Display the extracted surface")) +{ + addUpdateCallback("updateMesh", {&d_step, &d_IsoValue, &d_gridMin, &d_gridMax}, [this](const sofa::core::DataTracker&) + { + checkInputs(); + updateMeshIfNeeded(); + return core::objectmodel::ComponentState::Valid; + }, {&d_outPoints, &d_outTriangles}); + d_outPoints.setGroup("Output"); + d_outTriangles.setGroup("Output"); +} + +FieldToSurfaceMesh::~FieldToSurfaceMesh() +{ +} + +void FieldToSurfaceMesh::init() +{ + if(!l_field.get()) + { + msg_error() << "Missing field to extract surface from"; + d_componentState = core::objectmodel::ComponentState::Invalid; + } + + updateMeshIfNeeded(); + d_componentState = core::objectmodel::ComponentState::Valid; +} + +void FieldToSurfaceMesh::computeBBox(const core::ExecParams* /* params */, bool /*onlyVisible*/) +{ + f_bbox.setValue({d_gridMin.getValue(), d_gridMax.getValue()}); +} + +void FieldToSurfaceMesh::checkInputs(){ + + auto length = d_gridMax.getValue()-d_gridMin.getValue() ; + auto step = d_step.getValue(); + + // clamp the mStep value to avoid too large grids + if( step < 0.0001 || (length.x() / step > 256) || length.y() / step > 256 || length.z() / step > 256) + { + d_step.setValue( *std::max_element(length.begin(), length.end()) / 256.0 ); + msg_warning() << "step exceeding grid size, clamped to " << d_step.getValue(); + } +} + +void FieldToSurfaceMesh::updateMeshIfNeeded() +{ + sofa::helper::getWriteOnlyAccessor(d_outPoints).clear(); + sofa::helper::getWriteOnlyAccessor(d_outTriangles).clear(); + + double isoval = d_IsoValue.getValue(); + double mstep = d_step.getValue(); + double invStep = 1.0/d_step.getValue(); + + Vec3d gridmin = d_gridMin.getValue() ; + Vec3d gridmax = d_gridMax.getValue() ; + + auto field = l_field.get(); + + if(!field) + return; + + // Clear the previously used buffer + tmpPoints.clear(); + tmpTriangles.clear(); + + marchingCube.generateSurfaceMesh(isoval, mstep, invStep, gridmin, gridmax, + [field](const std::vector& positions, std::vector& res){ + field->getValues(positions, res); + }, + tmpPoints, tmpTriangles); + + /// Copy the surface to Sofa topology + d_outPoints.setValue(tmpPoints); + d_outTriangles.setValue(tmpTriangles); + + tmpPoints.clear(); + tmpTriangles.clear(); + + hasChanged = false; + return; +} + +void FieldToSurfaceMesh::draw(const VisualParams* vparams) +{ + if(isComponentStateInvalid()) + return; + + if(!d_debugDraw.getValue()) + return; + + auto drawTool = vparams->drawTool(); + + sofa::helper::ReadAccessor< Data > x = d_outPoints; + sofa::helper::ReadAccessor< Data > triangles = d_outTriangles; + drawTool->setLightingEnabled(true); + + for(const Triangle& triangle : triangles) + { + int a = triangle[0]; + int b = triangle[1]; + int c = triangle[2]; + Vec3d center = (x[a]+x[b]+x[c])*0.333333; + Vec3d pa = (0.9*x[a]+0.1*center) ; + Vec3d pb = (0.9*x[b]+0.1*center) ; + Vec3d pc = (0.9*x[c]+0.1*center) ; + + vparams->drawTool()->drawTriangles({pb,pa,pc}, + type::RGBAColor(0.0,0.0,1.0,1.0)); + } + + if(x.size()>1000){ + drawTool->drawPoints(x, 1.0, type::RGBAColor(1.0,1.0,0.0,0.2)) ; + }else{ + drawTool->drawSpheres(x, 0.01, type::RGBAColor(1.0,1.0,0.0,0.2)) ; + } +} + +// Register in the Factory +void registerFieldToSurfaceMesh(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Generates a surface mesh from a field function.") + .add< FieldToSurfaceMesh >()); +} + +} diff --git a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h new file mode 100644 index 00000000000..c4dbc73882e --- /dev/null +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h @@ -0,0 +1,102 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture, development version * +* (c) 2006-2025 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include +#include +#include +#include + +//////////////////////////////////////////////////////////////////////////////////////////////////// +namespace sofaimplicitfield::component::engine +{ +using namespace sofa; + +typedef sofa::core::topology::BaseMeshTopology::SeqTriangles SeqTriangles; +typedef sofa::core::topology::BaseMeshTopology::Triangle Triangle; +typedef sofa::type::vector VecCoord; + +using sofa::component::geometry::ScalarField; +using sofa::core::visual::VisualParams ; +using sofa::core::objectmodel::BaseObject ; +using sofa::type::Vec3d ; + +class FieldToSurfaceMesh : public BaseObject +{ +public: + SOFA_CLASS(FieldToSurfaceMesh, BaseObject); + + virtual void init() override ; + virtual void draw(const VisualParams*params) override ; + + double getStep() const { return d_step.getValue(); } + void setStep(double val) { d_step.setValue(val); } + + double getIsoValue() const { return d_IsoValue.getValue(); } + void setIsoValue(double val) { d_IsoValue.setValue(val); } + + const Vec3d& getGridMin() const { return d_gridMin.getValue(); } + void setGridMin(const Vec3d& val) { d_gridMin.setValue(val); } + void setGridMin(double x, double y, double z) { d_gridMin.setValue( Vec3d(x,y,z)); } + + const Vec3d& getGridMax() const { return d_gridMax.getValue(); } + void setGridMax(const Vec3d& val) { d_gridMax.setValue(val); } + void setGridMax(double x, double y, double z) { d_gridMax.setValue( Vec3d(x,y,z)); } + +protected: + SingleLink l_field ; + + Data d_step; + Data d_IsoValue; + + Data< Vec3d > d_gridMin; + Data< Vec3d > d_gridMax; + + /// Output + Data d_outPoints; + Data d_outTriangles; + Data d_debugDraw; + +protected: + FieldToSurfaceMesh() ; + virtual ~FieldToSurfaceMesh() ; + +private: + void computeBBox(const core::ExecParams* /* params */, bool /*onlyVisible*/=false) override; + + void checkInputs(); + + void generateSurfaceMesh(double isoval, double mstep, double invStep, + Vec3d gridmin, Vec3d gridmax, + sofa::component::geometry::ScalarField*); + void updateMeshIfNeeded(); + + bool hasChanged {true} ; + VecCoord tmpPoints; + SeqTriangles tmpTriangles; + + MarchingCube marchingCube; +}; + +} + diff --git a/applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp b/applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp index 178d034a6e6..2eae46b226b 100644 --- a/applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp +++ b/applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp @@ -40,6 +40,22 @@ namespace geometry namespace _scalarfield_ { +void ScalarField::init() +{ + d_componentState.setValue(core::objectmodel::ComponentState::Valid); +} + +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)); + } + return; +} + Vec3d ScalarField::getGradientByFinitDifference(Vec3d& pos, int& i) { Vec3d Result; diff --git a/applications/plugins/SofaImplicitField/components/geometry/ScalarField.h b/applications/plugins/SofaImplicitField/components/geometry/ScalarField.h index 51f28050e08..aa27b5dca04 100644 --- a/applications/plugins/SofaImplicitField/components/geometry/ScalarField.h +++ b/applications/plugins/SofaImplicitField/components/geometry/ScalarField.h @@ -49,6 +49,8 @@ class SOFA_SOFAIMPLICITFIELD_API ScalarField : public BaseObject SOFA_CLASS(ScalarField, BaseObject); public: + void init() override; + /// Compute the gradient using a first order finite-difference scheme. /// This is of lower precision compared to analytical gradient computed by derivating /// the equations. @@ -64,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); diff --git a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.cpp b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.cpp index 0b0c661042b..9609e654e7f 100644 --- a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.cpp +++ b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.cpp @@ -24,7 +24,7 @@ #include #include "ImplicitSurfaceMapping.inl" -namespace sofa::component::mapping +namespace sofaimplicitfield::mapping { using namespace sofa::defaulttype; @@ -38,6 +38,4 @@ void registerImplicitSurfaceMapping(sofa::core::ObjectFactory* factory) template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< Vec3dTypes, Vec3dTypes >; - -} // namespace sofa::component::mapping - +} diff --git a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.h b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.h index d4425812ea0..6b66ece9bd8 100644 --- a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.h +++ b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.h @@ -19,30 +19,26 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ -#ifndef SOFA_COMPONENT_MAPPING_IMPLICITSURFACEMAPPING_H -#define SOFA_COMPONENT_MAPPING_IMPLICITSURFACEMAPPING_H +#pragma once #include #include #include #include +#include #include - -namespace sofa -{ - -namespace component +namespace sofaimplicitfield::mapping { -namespace mapping -{ +using namespace sofa; +using sofa::component::topology::container::constant::MeshTopology; template -class ImplicitSurfaceMapping : public core::Mapping, public topology::container::constant::MeshTopology +class ImplicitSurfaceMapping : public core::Mapping, public MeshTopology { public: - SOFA_CLASS2(SOFA_TEMPLATE2(ImplicitSurfaceMapping, In, Out), SOFA_TEMPLATE2(core::Mapping, In, Out), topology::container::constant::MeshTopology); + SOFA_CLASS2(SOFA_TEMPLATE2(ImplicitSurfaceMapping, In, Out), SOFA_TEMPLATE2(core::Mapping, In, Out), MeshTopology); typedef core::Mapping Inherit; typedef typename Out::VecCoord OutVecCoord; @@ -114,6 +110,8 @@ class ImplicitSurfaceMapping : public core::Mapping, public topology::c msg_error() << "applyJT(constraint) is not implemented"; } + void draw(const core::visual::VisualParams* params) override; + protected: Data mStep; ///< Step Data mRadius; ///< Radius @@ -122,6 +120,10 @@ class ImplicitSurfaceMapping : public core::Mapping, public topology::c Data< InCoord > mGridMin; ///< Grid Min Data< InCoord > mGridMax; ///< Grid Max + Vec3d mLocalGridMin; ///< Grid Min + Vec3d mLocalGridMax; ///< Grid Max + + // Marching cube data /// For each cube, store the vertex indices on each 3 first edges, and the data value @@ -146,85 +148,17 @@ class ImplicitSurfaceMapping : public core::Mapping, public topology::c Data < sofa::type::vector > planes; typename sofa::type::vector::iterator P0; /// Pointer to first plane typename sofa::type::vector::iterator P1; /// Pointer to second plane - - void newPlane(); - - template - int addPoint(OutVecCoord& out, int x,int y,int z, OutReal v0, OutReal v1, OutReal iso) - { - int p = out.size(); - OutCoord pos = OutCoord((OutReal)x,(OutReal)y,(OutReal)z); - pos[C] -= (iso-v0)/(v1-v0); - out.resize(p+1); - out[p] = pos * mStep.getValue(); - return p; - } - - int addFace(int p1, int p2, int p3, int nbp) - { - if ((unsigned)p1<(unsigned)nbp && - (unsigned)p2<(unsigned)nbp && - (unsigned)p3<(unsigned)nbp) - { - SeqTriangles& triangles = *d_seqTriangles.beginEdit(); - int f = triangles.size(); - triangles.push_back(Triangle(p1, p3, p2)); - d_seqTriangles.endEdit(); - return f; - } - else - { - msg_error() << "Invalid face "< X - 11 / 10 / - | 3 | 1 - |/ |/ - 3----2----2 - / - / -|_ -Y - -*/ - - #if !defined(SOFA_COMPONENT_MAPPING_IMPLICITSURFACEMAPPING_CPP) extern template class SOFA_SOFAIMPLICITFIELD_API ImplicitSurfaceMapping< defaulttype::Vec3dTypes, defaulttype::Vec3dTypes >; - - - #endif +} // namespace -} // namespace mapping - -} // namespace component - -} // namespace sofa - -#endif diff --git a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl index b09348c1963..267b73e3fc3 100644 --- a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl +++ b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl @@ -19,8 +19,7 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ -#ifndef SOFA_COMPONENT_MAPPING_IMPLICITSURFACEMAPPING_INL -#define SOFA_COMPONENT_MAPPING_IMPLICITSURFACEMAPPING_INL +#pragma once #include "ImplicitSurfaceMapping.h" #include @@ -28,22 +27,14 @@ #include #include - - -namespace sofa -{ - -namespace component -{ - -namespace mapping +namespace sofaimplicitfield::mapping { template void ImplicitSurfaceMapping::init() { core::Mapping::init(); - topology::container::constant::MeshTopology::init(); + MeshTopology::init(); } template @@ -66,218 +57,108 @@ Real sqr(Real r) return r*r; } +template +void ImplicitSurfaceMapping::draw(const core::visual::VisualParams* params) +{ + auto dt = params->drawTool(); + + dt->drawBoundingBox(mGridMin.getValue(), mGridMax.getValue()); + dt->drawBoundingBox(mLocalGridMin, mLocalGridMax); +} + template void ImplicitSurfaceMapping::apply(const core::MechanicalParams * /*mparams*/, Data& dOut, const Data& dIn) { - OutVecCoord &out = *dOut.beginEdit(); const InVecCoord& in = dIn.getValue(); - InReal invStep = (InReal)(1/mStep.getValue()); - out.resize(0); clear(); if (in.size()==0) { + OutVecCoord &out = *dOut.beginEdit(); dOut.endEdit(); return; } - InReal xmin, xmax; - InReal ymin, ymax; - xmin = xmax = in[0][0]*invStep; - ymin = ymax = in[0][1]*invStep; - const InReal r = (InReal)(getRadius() / mStep.getValue()); - std::map > sortParticles; + auto minGrid = mGridMin.getValue(); + auto maxGrid = mGridMax.getValue(); + + InReal invStep = (InReal)(1/mStep.getValue()); + const InReal r = getRadius(); + + std::unordered_map > sortParticles; for (unsigned int ip=0; ip (*mGridMax.beginEdit())[0] || - c0[1] < (*mGridMin.beginEdit())[1] || c0[1] > (*mGridMax.beginEdit())[1] || - c0[2] < (*mGridMin.beginEdit())[2] || c0[2] > (*mGridMax.beginEdit())[2]) + if (c0[0] < minGrid[0] || c0[0] > maxGrid[0] || + c0[1] < minGrid[1] || c0[1] > maxGrid[1] || + c0[2] < minGrid[2] || c0[2] > maxGrid[2]) continue; - InCoord c = c0 * invStep; - if (c[0] < xmin) - xmin = c[0]; - else if (c[0] > xmax) - xmax = c[0]; - if (c[1] < ymin) - ymin = c[1]; - else if (c[1] > ymax) - ymax = c[1]; - int z0 = helper::rceil(c[2]-r); - int z1 = helper::rfloor(c[2]+r); - for (int z = z0; z < z1; ++z) + + InCoord c = c0 ; + int z0 = helper::rfloor((c[2]-r)*invStep); + int z1 = helper::rceil((c[2]+r)*invStep); + for (int z = z0; z <= z1; ++z) sortParticles[z].push_back(c); } - const int z0 = sortParticles.begin()->first - 1; - const int nz = sortParticles.rbegin()->first - z0 + 2; - const int y0 = helper::rceil(ymin-r) - 1; - const int ny = helper::rfloor(ymax+r) - y0 + 2; - const int x0 = helper::rceil(xmin-r) - 1; - const int nx = helper::rfloor(xmax+r) - x0 + 2; - - (*planes.beginEdit()).resize(2*nx*ny); - P0 = (*planes.beginEdit()).begin()+0; - P1 = (*planes.beginEdit()).begin()+nx*ny; + OutReal r2 = (OutReal)sqr(r); - //////// MARCHING CUBE //////// + double rr = getRadius(); + type::BoundingBox box{}; + for(auto& particle : in) + { + box.include(particle); + } + box.include(box.minBBox()+Vec3d{-rr,-rr,-rr}); + box.include(box.maxBBox()+Vec3d{+rr,+rr,+rr}); - const OutReal isoval = (OutReal) getIsoValue(); + mLocalGridMin = box.minBBox(); + mLocalGridMax = box.maxBBox(); - const int dx = 1; - const int dy = nx; - //const int dz = nx*ny; + type::BoundingBox bigBox {mGridMin.getValue(), mGridMax.getValue()}; + box.intersection(bigBox); - int x,y,z,i,mk; - const int *tri; + auto fieldFunction = [&sortParticles, &r, &r2, &invStep]( + std::vector& pos, std::vector& res) -> void { - OutReal r2 = (OutReal)sqr(r); - // First plane is all zero - z = 0; - newPlane(); - for (z=1; z& particles = sortParticles[z0+z]; - for (typename std::list::const_iterator it = particles.begin(); it != particles.end(); ++it) + int i = 0; + for(auto& position : pos ) { - InCoord c = *it; - int cx0 = helper::rceil(c[0]-r); - int cx1 = helper::rfloor(c[0]+r); - int cy0 = helper::rceil(c[1]-r); - int cy1 = helper::rfloor(c[1]+r); - OutCoord dp2; - dp2[2] = (OutReal)sqr(z0+z-c[2]); - i = (cx0-x0)+(cy0-y0)*nx; - for (int y = cy0 ; y <= cy1 ; y++) - { - dp2[1] = (OutReal)sqr(y-c[1]); - int ix = i; - for (int x = cx0 ; x <= cx1 ; x++, ix++) - { - dp2[0] = (OutReal)sqr(x-c[0]); - OutReal d2 = dp2[0]+dp2[1]+dp2[2]; - if (d2 < r2) - { - // Soft object field function from the Wyvill brothers - // See http://astronomy.swin.edu.au/~pbourke/modelling/implicitsurf/ - d2 /= r2; - (P1+ix)->data += (1 + (-4*d2*d2*d2 + 17*d2*d2 - 22*d2)/9); - } + double sumd = 0.0; + for(auto& particle : (particlesIt->second)){ + position.z() = z; + double d2 = (position - particle).norm2(); + if(d2 < r2){ + d2 /= r2; + sumd += (1 + (-4*d2*d2*d2 + 17*d2*d2 - 22*d2)/9); } - i += nx; } + res[i++] = sumd; } + return; + }; - i=0; - int edgecube[12]; - const int edgepts[12] = {0,1,0,1,0,1,0,1,2,2,2,2}; - typename std::vector::iterator base = (*planes.beginEdit()).begin(); - int ip0 = P0-base; - int ip1 = P1-base; - edgecube[0] = (ip0 -dy); - edgecube[1] = (ip0 ); - edgecube[2] = (ip0 ); - edgecube[3] = (ip0-dx ); - edgecube[4] = (ip1 -dy); - edgecube[5] = (ip1 ); - edgecube[6] = (ip1 ); - edgecube[7] = (ip1-dx ); - edgecube[8] = (ip1-dx-dy); - edgecube[9] = (ip1-dy ); - edgecube[10] = (ip1 ); - edgecube[11] = (ip1-dx ); - - // First line is all zero - { - y=0; - x=0; - i+=nx; - } - for(y=1; ydata>isoval)^((P1+i-dx)->data>isoval)) - { - (P1+i)->p[0] = addPoint<0>(out, x0+x,y0+y,z0+z,(P1+i)->data,(P1+i-dx)->data,isoval); - } - if (((P1+i)->data>isoval)^((P1+i-dy)->data>isoval)) - { - (P1+i)->p[1] = addPoint<1>(out, x0+x,y0+y,z0+z,(P1+i)->data,(P1+i-dy)->data,isoval); - } - if (((P1+i)->data>isoval)^((P0+i)->data>isoval)) - { - (P1+i)->p[2] = addPoint<2>(out, x0+x,y0+y,z0+z,(P1+i)->data,(P0+i)->data,isoval); - } - - // All points should now be created + auto triangles = helper::getWriteOnlyAccessor(d_seqTriangles); + auto points = helper::getWriteOnlyAccessor(dOut); - if ((P0+i-dx-dy)->data > isoval) mk = 1; else mk=0; - if ((P0+i -dy)->data > isoval) mk|= 2; - if ((P0+i )->data > isoval) mk|= 4; - if ((P0+i-dx )->data > isoval) mk|= 8; - if ((P1+i-dx-dy)->data > isoval) mk|= 16; - if ((P1+i -dy)->data > isoval) mk|= 32; - if ((P1+i )->data > isoval) mk|= 64; - if ((P1+i-dx )->data > isoval) mk|= 128; + points.clear(); + triangles.clear(); + marchingCube.generateSurfaceMesh(mIsoValue.getValue(), mStep.getValue(), + invStep, box.minBBox(), box.maxBBox(), + fieldFunction, points.wref(), triangles.wref()); - tri=sofa::helper::MarchingCubeTriTable[mk]; - while (*tri>=0) - { - typename std::vector::iterator b = base+i; - if (addFace((b+edgecube[tri[0]])->p[edgepts[tri[0]]], - (b+edgecube[tri[1]])->p[edgepts[tri[1]]], - (b+edgecube[tri[2]])->p[edgepts[tri[2]]], out.size())<0) - { - msg_error() << " mk=0x"< - + - + diff --git a/applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py b/applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py new file mode 100644 index 00000000000..02d506e2066 --- /dev/null +++ b/applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py @@ -0,0 +1,53 @@ +import Sofa +from Sofa.Types import RGBAColor +from xshape.primitives import * +from xshape.transforms import * +from xshape.operators import * + +class DrawController(Sofa.Core.Controller): + def __init__(self, *args, **kwargs): + Sofa.Core.Controller.__init__(self, *args, **kwargs) + + def draw(self, visual_context): + dt = visual_context.getDrawTool() + dt.drawText([-1.0, 1.0, 0.5], 0.2, "Union(Sphere, Box)", RGBAColor(1.0,1.0,1.0,1.0)) + dt.drawText([ 1.0, 1.0, 0.5], 0.2, "Difference(Sphere, Box)", RGBAColor(1.0,1.0,1.0,1.0)) + +def createScene(root : Sofa.Core.Node): + """Creates two different mesh from two scalar field. + 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", name="SofaImplicitField") + + root.addObject(DrawController()) + + ########################### Fields ################## + root.addChild("Fields") + f1 = root.Fields.addObject( + Union(name="field1", + childA=Sphere(name="sphere", center=[0,0,0],radius=0.7), + childB=RoundedBox(center=[0.0,0.0,0.0],dimensions=[0.95,0.5,0.5], rounding_radius=0.1)) + ) + + f2 = root.Fields.addObject( + Difference(name="field2", + childB=Sphere(name="sphere", center=[2,0,0],radius=0.9), + childA=RoundedBox(center=[2.0,0.0,0.0],dimensions=[0.95,0.5,0.5], rounding_radius=0.1)) + ) + + ########################### Meshing ################## + root.addChild("Meshing") + m1 = root.Meshing.addObject("FieldToSurfaceMesh", name="polygonizer1", + field=f1.linkpath, min=[-1,-1,-1], max=[1,1,1], + step=0.1, debugDraw=True) + + m2 = root.Meshing.addObject("FieldToSurfaceMesh", name="polygonizer2", + field=f2.linkpath, min=[1,-1,-1], max=[3,1,1], + step=0.07) + + ########################### Fields ################## + root.addChild("Visual") + root.Visual.addObject("OglModel", name="renderer", + position=root.Meshing.polygonizer2.points.linkpath, + triangles=root.Meshing.polygonizer2.triangles.linkpath) diff --git a/applications/plugins/SofaImplicitField/examples/python/python-implicit-field-example.py b/applications/plugins/SofaImplicitField/examples/python/python-implicit-field-example.py new file mode 100644 index 00000000000..ef82a608a09 --- /dev/null +++ b/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)) diff --git a/applications/plugins/SofaImplicitField/examples/python/xshape/__init__.py b/applications/plugins/SofaImplicitField/examples/python/xshape/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/applications/plugins/SofaImplicitField/examples/python/xshape/operators.py b/applications/plugins/SofaImplicitField/examples/python/xshape/operators.py new file mode 100644 index 00000000000..4506c1342a3 --- /dev/null +++ b/applications/plugins/SofaImplicitField/examples/python/xshape/operators.py @@ -0,0 +1,35 @@ +from SofaImplicitField import ScalarField +import numpy + +class Union(ScalarField): + """Union of two scalar fields""" + def __init__(self, *args, **kwargs): + ScalarField.__init__(self, *args, **kwargs) + + self.childA = kwargs.get("childA", None) + self.childB = kwargs.get("childB", None) + + def getValue(self, position): + return min(self.childA.getValue(position), self.childB.getValue(position)) + +class Difference(ScalarField): + """Difference of two scalar fields""" + def __init__(self, *args, **kwargs): + ScalarField.__init__(self, *args, **kwargs) + + self.childA = kwargs.get("childA", None) + self.childB = kwargs.get("childB", None) + + def getValue(self, position): + return max(-self.childA.getValue(position), self.childB.getValue(position)) + +class Intersection(ScalarField): + """Intersection of two scalar fields""" + def __init__(self, *args, **kwargs): + ScalarField.__init__(self, *args, **kwargs) + + self.childA = kwargs.get("childA", None) + 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 diff --git a/applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py b/applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py new file mode 100644 index 00000000000..fe203d79654 --- /dev/null +++ b/applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py @@ -0,0 +1,43 @@ +"""Distance field function + +Sources: + https://iquilezles.org/articles/distfunctions/ +""" +from SofaImplicitField import ScalarField +import numpy +import numba + +class Sphere(ScalarField): + def __init__(self, *args, **kwargs): + ScalarField.__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, 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) + + 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("dimensions", type="Vec3d",value=kwargs.get("dimensions", [1.0,1.0,1.0]), default=[1.0,1.0,1.0], help="dimmension of the box", group="Geometry") + self.addData("rounding_radius", type="double",value=kwargs.get("rounding_radius", 0.1), default=0.1, help="radius of the sphere", group="Geometry") + + def getValue(self, pos): + x,y,z = pos + b = self.dimensions.value + r = self.rounding_radius.value + 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 diff --git a/applications/plugins/SofaImplicitField/examples/python/xshape/transforms.py b/applications/plugins/SofaImplicitField/examples/python/xshape/transforms.py new file mode 100644 index 00000000000..778b4f7dd7c --- /dev/null +++ b/applications/plugins/SofaImplicitField/examples/python/xshape/transforms.py @@ -0,0 +1,15 @@ +from SofaImplicitField import ScalarField +import numpy + +class Translate(ScalarField): + """Translate a scalar field given as attribute""" + def __init__(self, *args, **kwargs): + ScalarField.__init__(self, *args, **kwargs) + + self.addData("translate", type="Vec3d",value=kwargs.get("translate", [0.0,0.0,0.0]), default=[0.0,0.0,0.0], help="amount of translation", group="Geometry") + self.child = kwargs.get("child", None) + + def getValue(self, pos): + x,y,z = pos + position = numpy.array([x,y,z])-self.translate.value + return self.child.getValue(position) \ No newline at end of file diff --git a/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp b/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp index ca5f0f5732d..9e3a2d51e6f 100644 --- a/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp +++ b/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp @@ -38,7 +38,7 @@ namespace sofa::component::geometry::_StarShapedField_ { extern void registerStarShapedField(sofa::core::ObjectFactory* factory); } -namespace sofa::component::mapping +namespace sofaimplicitfield::mapping { extern void registerImplicitSurfaceMapping(sofa::core::ObjectFactory* factory); } @@ -50,7 +50,10 @@ namespace sofa::component::geometry::_discretegrid_ { extern void registerDiscreteGridField(sofa::core::ObjectFactory* factory); } - +namespace sofaimplicitfield::component::engine +{ +extern void registerFieldToSurfaceMesh(sofa::core::ObjectFactory* factory); +} namespace sofaimplicitfield { @@ -103,9 +106,10 @@ void registerObjects(sofa::core::ObjectFactory* factory) sofa::component::geometry::_BottleField_::registerBottleField(factory); sofa::component::geometry::_sphericalfield_::registerSphericalField(factory); sofa::component::geometry::_StarShapedField_::registerStarShapedField(factory); - sofa::component::mapping::registerImplicitSurfaceMapping(factory); + sofaimplicitfield::mapping::registerImplicitSurfaceMapping(factory); sofa::component::container::registerInterpolatedImplicitSurface(factory); sofa::component::geometry::_discretegrid_::registerDiscreteGridField(factory); + sofaimplicitfield::component::engine::registerFieldToSurfaceMesh(factory); } } /// sofaimplicitfield diff --git a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp index fc2138831e4..30b7510c750 100644 --- a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp +++ b/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 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 +110,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); @@ -112,6 +182,7 @@ void moduleAddScalarField(py::module &m) { "positional argument='" + py::cast(args[0]) + "'."); } + ff->setName(py::cast(value)); } } return ff;