From 8144c9e54369de0cf8f412f4be49d842e35c4d73 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 27 Aug 2025 14:39:08 +0200 Subject: [PATCH 01/22] First version. --- .../plugins/SofaImplicitField/CMakeLists.txt | 2 + .../components/engine/FieldToSurfaceMesh.cpp | 317 ++++++++++++++++++ .../components/engine/FieldToSurfaceMesh.h | 131 ++++++++ .../components/geometry/ScalarField.cpp | 5 + .../components/geometry/ScalarField.h | 2 + .../example-mesh-extraction-from-implicit.py | 30 ++ .../examples/python/primitives.py | 13 + .../python/python-implicit-field-example.py | 17 + .../examples/python/python-scalarfield.py | 28 -- .../initSofaImplicitField.cpp | 7 +- 10 files changed, 523 insertions(+), 29 deletions(-) create mode 100644 applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp create mode 100644 applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h create mode 100644 applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py create mode 100644 applications/plugins/SofaImplicitField/examples/python/primitives.py create mode 100644 applications/plugins/SofaImplicitField/examples/python/python-implicit-field-example.py delete mode 100644 applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py diff --git a/applications/plugins/SofaImplicitField/CMakeLists.txt b/applications/plugins/SofaImplicitField/CMakeLists.txt index 35d08e7f7b1..d1fc7604abd 100644 --- a/applications/plugins/SofaImplicitField/CMakeLists.txt +++ b/applications/plugins/SofaImplicitField/CMakeLists.txt @@ -12,6 +12,7 @@ set(HEADER_FILES 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 @@ -28,6 +29,7 @@ set(SOURCE_FILES 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/components/engine/FieldToSurfaceMesh.cpp b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp new file mode 100644 index 00000000000..80b50d99ad9 --- /dev/null +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp @@ -0,0 +1,317 @@ +/****************************************************************************** +* 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.")) + , mStep(initData(&mStep,0.1,"step","Step")) + , mIsoValue(initData(&mIsoValue,0.0,"isoValue","Iso Value")) + , mGridMin(initData(&mGridMin, Vec3d(-1,-1,-1),"min","Grid Min")) + , mGridMax(initData(&mGridMax, 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", {&mStep, &mIsoValue, &mGridMin, &mGridMax}, [this](const sofa::core::DataTracker&) + { + checkInputs(); + hasChanged=true; + return core::objectmodel::ComponentState::Valid; + }, {}); +} + +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::checkInputs(){ + + auto length = mGridMax.getValue()-mGridMin.getValue() ; + auto step = mStep.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) + { + mStep.setValue( *std::max_element(length.begin(), length.end()) / 256.0 ); + msg_warning() << "step exceeding grid size, clamped to " << mStep.getValue(); + } +} + +void FieldToSurfaceMesh::updateMeshIfNeeded() +{ + if(!hasChanged) + return; + + sofa::helper::getWriteOnlyAccessor(d_outPoints).clear(); + sofa::helper::getWriteOnlyAccessor(d_outTriangles).clear(); + + double isoval = mIsoValue.getValue(); + double mstep = mStep.getValue(); + double invStep = 1.0/mStep.getValue(); + + Vec3d gridmin = mGridMin.getValue() ; + Vec3d gridmax = mGridMax.getValue() ; + + auto field = l_field.get(); + + generateSurfaceMesh(isoval, mstep, invStep, gridmin, gridmax, field); + + /// 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; + + updateMeshIfNeeded(); + + auto drawTool = vparams->drawTool(); + + drawTool->drawBoundingBox(mGridMin.getValue(), mGridMax.getValue()) ; + + 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) ; + + Vec3d a1 = x[c]-x[b] ; + Vec3d a2 = x[a]-x[b] ; + + vparams->drawTool()->drawTriangles({pa,pb,pc}, + a1.cross(a2), + 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)) ; + } +} + +void FieldToSurfaceMesh::generateSurfaceMesh(double isoval, double mstep, double invStep, + Vec3d gridmin, Vec3d gridmax, + sofa::component::geometry::ScalarField* field) +{ + if(!field) + return; + + tmpPoints.clear(); + tmpTriangles.clear(); + + 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 ; + + double cx,cy,cz; + int x,y,z,i,mk; + const int *tri; + + + planes.resize(2*(nx)*(ny)); + P0 = planes.begin()+0; + P1 = planes.begin()+nx*ny; + + const int dx = 1; + const int dy = nx; + + z = 0; + newPlane(); + + i = 0 ; + cz = gridmin.z() ; + for (int y = 0 ; y < ny ; ++y) + { + cy = gridmin.y() + mstep * y ; + for (int x = 0 ; x < nx ; ++x, ++i) + { + cx = gridmin.x() + mstep * x ; + + Vec3d pos { cx, cy, cz } ; + double res = field->getValue(pos) ; + (P1+i)->data = res ; + } + } + + for (z=1; z<=nz; ++z) + { + newPlane(); + + i = 0 ; + cz = gridmin.z() + mstep * z ; + for (int y = 0 ; y < ny ; ++y) + { + cy = gridmin.y() + mstep * y ; + for (int x = 0 ; x < nx ; ++x, ++i) + { + cx = gridmin.x() + mstep * x ; + + Vec3d pos { cx, cy, cz } ; + double res = field->getValue(pos) ; + (P1+i)->data = res ; + } + } + + unsigned int 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.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(tmpPoints, 0, pos,gridmin, (P1+i)->data,(P1+i-dx)->data, mstep, isoval); + } + if (((P1+i)->data>isoval)^((P1+i-dy)->data>isoval)) + { + (P1+i)->p[1] = addPoint(tmpPoints, 1, pos,gridmin,(P1+i)->data,(P1+i-dy)->data, mstep, isoval); + } + if (((P1+i)->data>isoval)^((P0+i)->data>isoval)) + { + (P1+i)->p[2] = addPoint(tmpPoints, 2, pos,gridmin,(P1+i)->data,(P0+i)->data, mstep, isoval); + } + + // All points should now be created + 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; + + tri=sofa::helper::MarchingCubeTriTable[mk]; + while (*tri>=0) + { + typename std::vector::iterator b = base+i; + 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; + } + ++i; + } + } + } +} + +void FieldToSurfaceMesh::newPlane() +{ + CubeData c; + c.p[0] = -1; + c.p[1] = -1; + c.p[2] = -1; + c.data = 0; + typename std::vector::iterator P = P0; + P0 = P1; + P1 = P; + int n = planes.size()/2; + for (int i=0; iregisterObjects(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..aa49a52d8ab --- /dev/null +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h @@ -0,0 +1,131 @@ +/****************************************************************************** +* 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::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 mStep.getValue(); } + void setStep(double val) { mStep.setValue(val); } + + double getIsoValue() const { return mIsoValue.getValue(); } + void setIsoValue(double val) { mIsoValue.setValue(val); } + + const Vec3d& getGridMin() const { return mGridMin.getValue(); } + void setGridMin(const Vec3d& val) { mGridMin.setValue(val); } + void setGridMin(double x, double y, double z) { mGridMin.setValue( Vec3d(x,y,z)); } + + const Vec3d& getGridMax() const { return mGridMax.getValue(); } + void setGridMax(const Vec3d& val) { mGridMax.setValue(val); } + void setGridMax(double x, double y, double z) { mGridMax.setValue( Vec3d(x,y,z)); } + +protected: + SingleLink l_field ; + + Data mStep; + Data mIsoValue; + + Data< Vec3d > mGridMin; + Data< Vec3d > mGridMax; + + /// For each cube, store the vertex indices on each 3 first edges, and the data value + struct CubeData + { + int p[3]; + double data; + }; + + int addPoint(VecCoord& 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; + } + } + + /// Output + Data d_outPoints; + Data d_outTriangles; + Data d_debugDraw; + + sofa::type::vector planes; + typename sofa::type::vector::iterator P0; /// Pointer to first plane + typename sofa::type::vector::iterator P1; /// Pointer to second plane + +protected: + FieldToSurfaceMesh() ; + virtual ~FieldToSurfaceMesh() ; + +private: + void checkInputs(); + void newPlane(); + 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; +}; + +} + diff --git a/applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp b/applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp index 178d034a6e6..9c0323bbfd8 100644 --- a/applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp +++ b/applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp @@ -40,6 +40,11 @@ namespace geometry namespace _scalarfield_ { +void ScalarField::init() +{ + d_componentState.setValue(core::objectmodel::ComponentState::Valid); +} + 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..bebd099e6f0 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. 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..895a7191165 --- /dev/null +++ b/applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py @@ -0,0 +1,30 @@ +import Sofa +from primitives import Sphere + +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") + + ########################### Fields ################## + root.addChild("Fields") + f1 = root.Fields.addObject(Sphere(name="field1", center=[0,0,0])) + f2 = root.Fields.addObject("SphericalField", name="field2", center=[2,0,0]) + + ########################### 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.01) + + ########################### 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/primitives.py b/applications/plugins/SofaImplicitField/examples/python/primitives.py new file mode 100644 index 00000000000..b7b8fb96908 --- /dev/null +++ b/applications/plugins/SofaImplicitField/examples/python/primitives.py @@ -0,0 +1,13 @@ + +from SofaImplicitField import ScalarField +import numpy + +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, x, y, z): + return numpy.sqrt( numpy.sum((self.center.value - numpy.array([x,y,z]))**2) ) - self.radius.value \ No newline at end of file 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..6612ffb2fca --- /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/python-scalarfield.py b/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py deleted file mode 100644 index d97e36a5814..00000000000 --- a/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py +++ /dev/null @@ -1,28 +0,0 @@ -import Sofa -from SofaImplicitField import ScalarField -import numpy - -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, x, y, z): - return numpy.sqrt( numpy.sum((self.center.value - numpy.array([x,y,z]))**2) ) - self.radius.value - -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/initSofaImplicitField.cpp b/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp index ca5f0f5732d..53c301dbb02 100644 --- a/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp +++ b/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp @@ -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 { @@ -100,12 +103,14 @@ const char* getModuleDescription() void registerObjects(sofa::core::ObjectFactory* factory) { + std::cout << "===================================================================" << std::endl; sofa::component::geometry::_BottleField_::registerBottleField(factory); sofa::component::geometry::_sphericalfield_::registerSphericalField(factory); sofa::component::geometry::_StarShapedField_::registerStarShapedField(factory); sofa::component::mapping::registerImplicitSurfaceMapping(factory); sofa::component::container::registerInterpolatedImplicitSurface(factory); sofa::component::geometry::_discretegrid_::registerDiscreteGridField(factory); + sofaimplicitfield::component::engine::registerFieldToSurfaceMesh(factory); } } /// sofaimplicitfield From 3623617b6a00f8928d5e6bbac83b5f2106338903 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 3 Sep 2025 14:12:36 +0200 Subject: [PATCH 02/22] Add a MarchingCube implementation in a separated file. --- .../plugins/SofaImplicitField/CMakeLists.txt | 2 + .../SofaImplicitField/MarchingCube.cpp | 185 ++++++++++++++++++ .../plugins/SofaImplicitField/MarchingCube.h | 85 ++++++++ 3 files changed, 272 insertions(+) create mode 100644 applications/plugins/SofaImplicitField/MarchingCube.cpp create mode 100644 applications/plugins/SofaImplicitField/MarchingCube.h diff --git a/applications/plugins/SofaImplicitField/CMakeLists.txt b/applications/plugins/SofaImplicitField/CMakeLists.txt index d1fc7604abd..aca4a06b47b 100644 --- a/applications/plugins/SofaImplicitField/CMakeLists.txt +++ b/applications/plugins/SofaImplicitField/CMakeLists.txt @@ -6,6 +6,7 @@ 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 @@ -24,6 +25,7 @@ set(HEADER_FILES set(SOURCE_FILES initSofaImplicitField.cpp + MarchingCube.cpp ## This is a backward compatibility.. deprecated/SphereSurface.cpp diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp new file mode 100644 index 00000000000..774125b2d84 --- /dev/null +++ b/applications/plugins/SofaImplicitField/MarchingCube.cpp @@ -0,0 +1,185 @@ +/****************************************************************************** +* 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::newPlane() +{ + CubeData c; + c.p[0] = -1; + c.p[1] = -1; + c.p[2] = -1; + c.data = 0; + typename std::vector::iterator P = P0; + P0 = P1; + P1 = P; + int n = planes.size()/2; + for (int i=0; igetValue(pos) ; + (P1+i)->data = res ; + } + } + + for (z=1; z<=nz; ++z) + { + newPlane(); + + i = 0 ; + cz = gridmin.z() + mstep * z ; + for (int y = 0 ; y < ny ; ++y) + { + cy = gridmin.y() + mstep * y ; + for (int x = 0 ; x < nx ; ++x, ++i) + { + cx = gridmin.x() + mstep * x ; + + Vec3d pos { cx, cy, cz } ; + double res = field->getValue(pos) ; + (P1+i)->data = res ; + } + } + + unsigned int 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.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(tmpPoints, 0, pos,gridmin, (P1+i)->data,(P1+i-dx)->data, mstep, isoval); + } + if (((P1+i)->data>isoval)^((P1+i-dy)->data>isoval)) + { + (P1+i)->p[1] = addPoint(tmpPoints, 1, pos,gridmin,(P1+i)->data,(P1+i-dy)->data, mstep, isoval); + } + if (((P1+i)->data>isoval)^((P0+i)->data>isoval)) + { + (P1+i)->p[2] = addPoint(tmpPoints, 2, pos,gridmin,(P1+i)->data,(P0+i)->data, mstep, isoval); + } + + // All points should now be created + 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; + + tri=sofa::helper::MarchingCubeTriTable[mk]; + while (*tri>=0) + { + typename std::vector::iterator b = base+i; + 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; + } + ++i; + } + } + } +} + +} diff --git a/applications/plugins/SofaImplicitField/MarchingCube.h b/applications/plugins/SofaImplicitField/MarchingCube.h new file mode 100644 index 00000000000..91fac485258 --- /dev/null +++ b/applications/plugins/SofaImplicitField/MarchingCube.h @@ -0,0 +1,85 @@ +/****************************************************************************** +* 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, + sofa::component::geometry::ScalarField* field, + SeqCoord& tmpPoints, SeqTriangles& tmpTriangles); + +private: + void newPlane(); + + /// 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; + } + } +}; + + +} + From c31b1cd5d612562f1d3172fedcddcba23993f39b Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 3 Sep 2025 14:13:53 +0200 Subject: [PATCH 03/22] Use MarchingCube in FieldToSurfaceMesh --- .../components/engine/FieldToSurfaceMesh.cpp | 191 +++--------------- .../components/engine/FieldToSurfaceMesh.h | 69 ++----- 2 files changed, 42 insertions(+), 218 deletions(-) diff --git a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp index 80b50d99ad9..b486c9138f4 100644 --- a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp @@ -36,15 +36,15 @@ namespace sofaimplicitfield::component::engine FieldToSurfaceMesh::FieldToSurfaceMesh() : l_field(initLink("field", "The scalar field to generate a mesh from.")) - , mStep(initData(&mStep,0.1,"step","Step")) - , mIsoValue(initData(&mIsoValue,0.0,"isoValue","Iso Value")) - , mGridMin(initData(&mGridMin, Vec3d(-1,-1,-1),"min","Grid Min")) - , mGridMax(initData(&mGridMax, Vec3d(1,1,1),"max","Grid Max")) + , 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_dridMax(initData(&d_dridMax, 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", {&mStep, &mIsoValue, &mGridMin, &mGridMax}, [this](const sofa::core::DataTracker&) + addUpdateCallback("updateMesh", {&d_step, &d_IsoValue, &d_gridMin, &d_dridMax}, [this](const sofa::core::DataTracker&) { checkInputs(); hasChanged=true; @@ -71,14 +71,14 @@ void FieldToSurfaceMesh::init() void FieldToSurfaceMesh::checkInputs(){ - auto length = mGridMax.getValue()-mGridMin.getValue() ; - auto step = mStep.getValue(); + auto length = d_dridMax.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) { - mStep.setValue( *std::max_element(length.begin(), length.end()) / 256.0 ); - msg_warning() << "step exceeding grid size, clamped to " << mStep.getValue(); + d_step.setValue( *std::max_element(length.begin(), length.end()) / 256.0 ); + msg_warning() << "step exceeding grid size, clamped to " << d_step.getValue(); } } @@ -90,16 +90,23 @@ void FieldToSurfaceMesh::updateMeshIfNeeded() sofa::helper::getWriteOnlyAccessor(d_outPoints).clear(); sofa::helper::getWriteOnlyAccessor(d_outTriangles).clear(); - double isoval = mIsoValue.getValue(); - double mstep = mStep.getValue(); - double invStep = 1.0/mStep.getValue(); + double isoval = d_IsoValue.getValue(); + double mstep = d_step.getValue(); + double invStep = 1.0/d_step.getValue(); - Vec3d gridmin = mGridMin.getValue() ; - Vec3d gridmax = mGridMax.getValue() ; + Vec3d gridmin = d_gridMin.getValue() ; + Vec3d gridmax = d_dridMax.getValue() ; auto field = l_field.get(); - generateSurfaceMesh(isoval, mstep, invStep, gridmin, gridmax, field); + if(!field) + return; + + // Clear the previously used buffer + tmpPoints.clear(); + tmpTriangles.clear(); + + marchingCube.generateSurfaceMesh(isoval, mstep, invStep, gridmin, gridmax, field, tmpPoints, tmpTriangles); /// Copy the surface to Sofa topology d_outPoints.setValue(tmpPoints); @@ -124,7 +131,7 @@ void FieldToSurfaceMesh::draw(const VisualParams* vparams) auto drawTool = vparams->drawTool(); - drawTool->drawBoundingBox(mGridMin.getValue(), mGridMax.getValue()) ; + drawTool->drawBoundingBox(d_gridMin.getValue(), d_dridMax.getValue()) ; sofa::helper::ReadAccessor< Data > x = d_outPoints; sofa::helper::ReadAccessor< Data > triangles = d_outTriangles; @@ -155,158 +162,6 @@ void FieldToSurfaceMesh::draw(const VisualParams* vparams) } } -void FieldToSurfaceMesh::generateSurfaceMesh(double isoval, double mstep, double invStep, - Vec3d gridmin, Vec3d gridmax, - sofa::component::geometry::ScalarField* field) -{ - if(!field) - return; - - tmpPoints.clear(); - tmpTriangles.clear(); - - 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 ; - - double cx,cy,cz; - int x,y,z,i,mk; - const int *tri; - - - planes.resize(2*(nx)*(ny)); - P0 = planes.begin()+0; - P1 = planes.begin()+nx*ny; - - const int dx = 1; - const int dy = nx; - - z = 0; - newPlane(); - - i = 0 ; - cz = gridmin.z() ; - for (int y = 0 ; y < ny ; ++y) - { - cy = gridmin.y() + mstep * y ; - for (int x = 0 ; x < nx ; ++x, ++i) - { - cx = gridmin.x() + mstep * x ; - - Vec3d pos { cx, cy, cz } ; - double res = field->getValue(pos) ; - (P1+i)->data = res ; - } - } - - for (z=1; z<=nz; ++z) - { - newPlane(); - - i = 0 ; - cz = gridmin.z() + mstep * z ; - for (int y = 0 ; y < ny ; ++y) - { - cy = gridmin.y() + mstep * y ; - for (int x = 0 ; x < nx ; ++x, ++i) - { - cx = gridmin.x() + mstep * x ; - - Vec3d pos { cx, cy, cz } ; - double res = field->getValue(pos) ; - (P1+i)->data = res ; - } - } - - unsigned int 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.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(tmpPoints, 0, pos,gridmin, (P1+i)->data,(P1+i-dx)->data, mstep, isoval); - } - if (((P1+i)->data>isoval)^((P1+i-dy)->data>isoval)) - { - (P1+i)->p[1] = addPoint(tmpPoints, 1, pos,gridmin,(P1+i)->data,(P1+i-dy)->data, mstep, isoval); - } - if (((P1+i)->data>isoval)^((P0+i)->data>isoval)) - { - (P1+i)->p[2] = addPoint(tmpPoints, 2, pos,gridmin,(P1+i)->data,(P0+i)->data, mstep, isoval); - } - - // All points should now be created - 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; - - tri=sofa::helper::MarchingCubeTriTable[mk]; - while (*tri>=0) - { - typename std::vector::iterator b = base+i; - 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; - } - ++i; - } - } - } -} - -void FieldToSurfaceMesh::newPlane() -{ - CubeData c; - c.p[0] = -1; - c.p[1] = -1; - c.p[2] = -1; - c.data = 0; - typename std::vector::iterator P = P0; - P0 = P1; - P1 = P; - int n = planes.size()/2; - for (int i=0; i #include - +#include #include #include @@ -48,75 +48,42 @@ class FieldToSurfaceMesh : public BaseObject virtual void init() override ; virtual void draw(const VisualParams*params) override ; - double getStep() const { return mStep.getValue(); } - void setStep(double val) { mStep.setValue(val); } + double getStep() const { return d_step.getValue(); } + void setStep(double val) { d_step.setValue(val); } - double getIsoValue() const { return mIsoValue.getValue(); } - void setIsoValue(double val) { mIsoValue.setValue(val); } + double getIsoValue() const { return d_IsoValue.getValue(); } + void setIsoValue(double val) { d_IsoValue.setValue(val); } - const Vec3d& getGridMin() const { return mGridMin.getValue(); } - void setGridMin(const Vec3d& val) { mGridMin.setValue(val); } - void setGridMin(double x, double y, double z) { mGridMin.setValue( Vec3d(x,y,z)); } + 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 mGridMax.getValue(); } - void setGridMax(const Vec3d& val) { mGridMax.setValue(val); } - void setGridMax(double x, double y, double z) { mGridMax.setValue( Vec3d(x,y,z)); } + const Vec3d& getGridMax() const { return d_dridMax.getValue(); } + void setGridMax(const Vec3d& val) { d_dridMax.setValue(val); } + void setGridMax(double x, double y, double z) { d_dridMax.setValue( Vec3d(x,y,z)); } protected: SingleLink l_field ; - Data mStep; - Data mIsoValue; - - Data< Vec3d > mGridMin; - Data< Vec3d > mGridMax; - - /// For each cube, store the vertex indices on each 3 first edges, and the data value - struct CubeData - { - int p[3]; - double data; - }; - - int addPoint(VecCoord& 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; - } - } + Data d_step; + Data d_IsoValue; + + Data< Vec3d > d_gridMin; + Data< Vec3d > d_dridMax; /// Output Data d_outPoints; Data d_outTriangles; Data d_debugDraw; - sofa::type::vector planes; - typename sofa::type::vector::iterator P0; /// Pointer to first plane - typename sofa::type::vector::iterator P1; /// Pointer to second plane - protected: FieldToSurfaceMesh() ; virtual ~FieldToSurfaceMesh() ; private: void checkInputs(); - void newPlane(); + void generateSurfaceMesh(double isoval, double mstep, double invStep, Vec3d gridmin, Vec3d gridmax, sofa::component::geometry::ScalarField*); @@ -125,6 +92,8 @@ class FieldToSurfaceMesh : public BaseObject bool hasChanged {true} ; VecCoord tmpPoints; SeqTriangles tmpTriangles; + + MarchingCube marchingCube; }; } From 6accc38e2eec28a4938555076dd1f48fcb450e1c Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 3 Sep 2025 14:15:59 +0200 Subject: [PATCH 04/22] Small fix-up on the signature for getValue in BInding_ScalarField --- .../plugins/SofaImplicitField/examples/python/primitives.py | 3 ++- .../SofaImplicitField/python/src/Binding_ScalarField.cpp | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/applications/plugins/SofaImplicitField/examples/python/primitives.py b/applications/plugins/SofaImplicitField/examples/python/primitives.py index b7b8fb96908..174e9ddbd7d 100644 --- a/applications/plugins/SofaImplicitField/examples/python/primitives.py +++ b/applications/plugins/SofaImplicitField/examples/python/primitives.py @@ -9,5 +9,6 @@ def __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, x, y, z): + def getValue(self, pos): + x,y,z = pos return numpy.sqrt( numpy.sum((self.center.value - numpy.array([x,y,z]))**2) ) - self.radius.value \ No newline at end of file diff --git a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp index 9e2243bd5ab..4110addc16e 100644 --- a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp +++ b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp @@ -43,7 +43,7 @@ class ScalarField_Trampoline : public ScalarField { SOFA_UNUSED(domain); PythonEnvironment::gil acquire; - PYBIND11_OVERLOAD_PURE(double, ScalarField, getValue, pos.x(), pos.y(), pos.z()); + PYBIND11_OVERLOAD_PURE(double, ScalarField, getValue, pos); } }; From 271999bf242f1a0ed237cde74b689eac3a49edf1 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 3 Sep 2025 14:44:02 +0200 Subject: [PATCH 05/22] Update to sync with the small fix --- .../examples/python/python-implicit-field-example.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/applications/plugins/SofaImplicitField/examples/python/python-implicit-field-example.py b/applications/plugins/SofaImplicitField/examples/python/python-implicit-field-example.py index 6612ffb2fca..ef82a608a09 100644 --- a/applications/plugins/SofaImplicitField/examples/python/python-implicit-field-example.py +++ b/applications/plugins/SofaImplicitField/examples/python/python-implicit-field-example.py @@ -8,9 +8,9 @@ def __init__(self, *args, **kwargs): 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) ) + 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")) From 7eeb3c811aae4320ef1055f745119281fa4464b7 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 3 Sep 2025 14:58:15 +0200 Subject: [PATCH 06/22] Use std::function in MarchingCube::generateField --- .../plugins/SofaImplicitField/MarchingCube.cpp | 12 +++--------- .../plugins/SofaImplicitField/MarchingCube.h | 2 +- .../components/engine/FieldToSurfaceMesh.cpp | 4 +++- 3 files changed, 7 insertions(+), 11 deletions(-) diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp index 774125b2d84..771de74f5f2 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.cpp +++ b/applications/plugins/SofaImplicitField/MarchingCube.cpp @@ -47,15 +47,9 @@ void MarchingCube::newPlane() void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, const double invStep, const Vec3d& gridmin, const Vec3d& gridmax, - sofa::component::geometry::ScalarField* field, + std::function getFieldValueAt, SeqCoord& tmpPoints, SeqTriangles& tmpTriangles) { - if(!field) - return; - - tmpPoints.clear(); - tmpTriangles.clear(); - 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 ; @@ -84,7 +78,7 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, cx = gridmin.x() + mstep * x ; Vec3d pos { cx, cy, cz } ; - double res = field->getValue(pos) ; + double res = getFieldValueAt(pos) ; (P1+i)->data = res ; } } @@ -103,7 +97,7 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, cx = gridmin.x() + mstep * x ; Vec3d pos { cx, cy, cz } ; - double res = field->getValue(pos) ; + double res = getFieldValueAt(pos) ; (P1+i)->data = res ; } } diff --git a/applications/plugins/SofaImplicitField/MarchingCube.h b/applications/plugins/SofaImplicitField/MarchingCube.h index 91fac485258..85e2e6d49f6 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.h +++ b/applications/plugins/SofaImplicitField/MarchingCube.h @@ -40,7 +40,7 @@ class MarchingCube public: void generateSurfaceMesh(const double isoval, const double mstep, const double invStep, const Vec3d& gridmin, const Vec3d& gridmax, - sofa::component::geometry::ScalarField* field, + std::function field, SeqCoord& tmpPoints, SeqTriangles& tmpTriangles); private: diff --git a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp index b486c9138f4..311e7816dd3 100644 --- a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp @@ -106,7 +106,9 @@ void FieldToSurfaceMesh::updateMeshIfNeeded() tmpPoints.clear(); tmpTriangles.clear(); - marchingCube.generateSurfaceMesh(isoval, mstep, invStep, gridmin, gridmax, field, tmpPoints, tmpTriangles); + marchingCube.generateSurfaceMesh(isoval, mstep, invStep, gridmin, gridmax, + [field](Vec3d& pos){return field->getValue(pos);}, + tmpPoints, tmpTriangles); /// Copy the surface to Sofa topology d_outPoints.setValue(tmpPoints); From f38934372923383960d15a5752bb88d7ca9fe7cc Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 3 Sep 2025 15:27:32 +0200 Subject: [PATCH 07/22] Refactoring ImplicitSurfaceMapping [WIP] --- .../SofaImplicitField/MarchingCube.cpp | 3 + .../mapping/ImplicitSurfaceMapping.cpp | 6 +- .../mapping/ImplicitSurfaceMapping.h | 92 +------- .../mapping/ImplicitSurfaceMapping.inl | 206 +++++------------- .../initSofaImplicitField.cpp | 4 +- 5 files changed, 72 insertions(+), 239 deletions(-) diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp index 771de74f5f2..5d6aae27c19 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.cpp +++ b/applications/plugins/SofaImplicitField/MarchingCube.cpp @@ -83,8 +83,11 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, } } + std::cout << "FIRST PLANE DONE " << std::endl; + for (z=1; z<=nz; ++z) { + std::cout << " PLANE DONE " << z << std::endl; newPlane(); i = 0 ; 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..d92e7a0a82d 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; @@ -146,85 +142,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..49fd084d785 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 @@ -135,149 +126,62 @@ void ImplicitSurfaceMapping::apply(const core::MechanicalParams * /*mpar 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) - { - 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); - } - } - i += nx; - } - } - - 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 - - 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; - - 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"<::apply(const core::MechanicalParams * /*mpar 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; - - //////// MARCHING CUBE //////// - - const OutReal isoval = (OutReal) getIsoValue(); - - const int dx = 1; - const int dy = nx; - //const int dz = nx*ny; - - int x,y,z,i,mk; - const int *tri; - 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) -// { -// 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); -// } -// } -// i += nx; -// } -// } -// } - - auto fieldFunction = [](Vec3d& pos) -> double { - return 0.5; + type::BoundingBox bigBox {mGridMin.getValue(), mGridMax.getValue()}; + type::BoundingBox box; + for(auto& [_, z_plane] : sortParticles) + { + for(auto& particle : z_plane) + box.include(particle); + } + box.intersection(bigBox); + + auto fieldFunction = [&sortParticles, &r, &r2](Vec3d& pos) -> double { + int index = helper::rceil(pos.z()); + + double sumd = 0.0; + for(auto& particle : sortParticles[index]){ + double d2 = (pos - particle).norm2(); + if(d2 < r2){ + d2 /= r2; + sumd += (1 + (-4*d2*d2*d2 + 17*d2*d2 - 22*d2)/9); + } + } + return sumd; }; - SeqTriangles triangles = helper::getWriteOnlyAccessor(d_seqTriangles); - SeqPoints points = helper::getWriteOnlyAccessor(dOut); + auto triangles = helper::getWriteOnlyAccessor(d_seqTriangles); + auto points = helper::getWriteOnlyAccessor(dOut); points.clear(); triangles.clear(); - marchingCube.generateSurfaceMesh(mIsoValue.getValue(), mStep.getValue(), - invStep, mGridMin.getValue(), mGridMax.getValue(), - fieldFunction, points, triangles); + invStep, box.minBBox(), box.maxBBox(), + fieldFunction, points.wref(), triangles.wref()); + + } template From f1cbc17c9a8842c53cc23daefc1a7358e0595ec7 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 3 Sep 2025 23:48:44 +0200 Subject: [PATCH 09/22] Improved refactoring of ImplicitSurfaceMapping --- .../SofaImplicitField/MarchingCube.cpp | 3 + .../mapping/ImplicitSurfaceMapping.h | 6 ++ .../mapping/ImplicitSurfaceMapping.inl | 75 +++++++++++-------- .../examples/ImplicitSurfaceMapping.scn | 7 +- 4 files changed, 57 insertions(+), 34 deletions(-) diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp index 771de74f5f2..bea66804a58 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.cpp +++ b/applications/plugins/SofaImplicitField/MarchingCube.cpp @@ -54,6 +54,9 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, int ny = floor((gridmax.y() - gridmin.y()) * invStep) + 1 ; int nz = floor((gridmax.z() - gridmin.z()) * invStep) + 1 ; + if( nz <= 0 || ny <= 0 || nx <= 0 ) + return; + double cx,cy,cz; int x,y,z,i,mk; const int *tri; diff --git a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.h b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.h index d92e7a0a82d..6b66ece9bd8 100644 --- a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.h +++ b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.h @@ -110,6 +110,8 @@ class ImplicitSurfaceMapping : public core::Mapping, public MeshTopolog msg_error() << "applyJT(constraint) is not implemented"; } + void draw(const core::visual::VisualParams* params) override; + protected: Data mStep; ///< Step Data mRadius; ///< Radius @@ -118,6 +120,10 @@ class ImplicitSurfaceMapping : public core::Mapping, public MeshTopolog 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 diff --git a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl index f15056c9fe2..f6467f631ef 100644 --- a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl +++ b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl @@ -57,65 +57,76 @@ 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); } OutReal r2 = (OutReal)sqr(r); - type::BoundingBox bigBox {mGridMin.getValue(), mGridMax.getValue()}; - type::BoundingBox box; - for(auto& [_, z_plane] : sortParticles) + + double rr = getRadius(); + type::BoundingBox box{}; + for(auto& particle : in) { - for(auto& particle : z_plane) - box.include(particle); + box.include(particle); } + box.include(box.minBBox()+Vec3d{-rr,-rr,-rr}); + box.include(box.maxBBox()+Vec3d{+rr,+rr,+rr}); + + mLocalGridMin = box.minBBox(); + mLocalGridMax = box.maxBBox(); + + type::BoundingBox bigBox {mGridMin.getValue(), mGridMax.getValue()}; box.intersection(bigBox); - auto fieldFunction = [&sortParticles, &r, &r2](Vec3d& pos) -> double { - int index = helper::rceil(pos.z()); + auto fieldFunction = [&sortParticles, &r, &r2, &invStep](Vec3d& pos) -> double { + int index = helper::rfloor(pos.z()*invStep); + auto particlesIt = sortParticles.find(index); + if(particlesIt==sortParticles.end()) + return 0.0; double sumd = 0.0; - for(auto& particle : sortParticles[index]){ + for(auto& particle : (particlesIt->second)){ double d2 = (pos - particle).norm2(); if(d2 < r2){ d2 /= r2; diff --git a/applications/plugins/SofaImplicitField/examples/ImplicitSurfaceMapping.scn b/applications/plugins/SofaImplicitField/examples/ImplicitSurfaceMapping.scn index 39da4815826..065f5ea055d 100644 --- a/applications/plugins/SofaImplicitField/examples/ImplicitSurfaceMapping.scn +++ b/applications/plugins/SofaImplicitField/examples/ImplicitSurfaceMapping.scn @@ -12,11 +12,14 @@ - + - + From e804f4a98fb1b4a3ee6e0c876bbf0b85e7f6b91c Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Thu, 4 Sep 2025 21:50:47 +0200 Subject: [PATCH 10/22] Make a "vectorized" version of the marching cube... --- .../SofaImplicitField/MarchingCube.cpp | 51 +++++++++++++------ .../plugins/SofaImplicitField/MarchingCube.h | 2 +- .../components/engine/FieldToSurfaceMesh.cpp | 8 ++- .../mapping/ImplicitSurfaceMapping.inl | 32 +++++++----- 4 files changed, 63 insertions(+), 30 deletions(-) diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp index bea66804a58..283ab21d3dc 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.cpp +++ b/applications/plugins/SofaImplicitField/MarchingCube.cpp @@ -47,7 +47,7 @@ void MarchingCube::newPlane() void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, const double invStep, const Vec3d& gridmin, const Vec3d& gridmax, - std::function getFieldValueAt, + std::function&, std::vector&)> getFieldValueAt, SeqCoord& tmpPoints, SeqTriangles& tmpTriangles) { int nx = floor((gridmax.x() - gridmin.x()) * invStep) + 1 ; @@ -73,36 +73,55 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, i = 0 ; cz = gridmin.z() ; + + std::vector positions; + std::vector output; + positions.resize(ny*nx); + output.resize(nx*ny); for (int y = 0 ; y < ny ; ++y) { cy = gridmin.y() + mstep * y ; for (int x = 0 ; x < nx ; ++x, ++i) { cx = gridmin.x() + mstep * x ; - - Vec3d pos { cx, cy, cz } ; - double res = getFieldValueAt(pos) ; - (P1+i)->data = res ; + positions[i].set(cx, cy, cz ); } } + getFieldValueAt(positions, output) ; + + // Copy back the data into planes. + auto it = P1; + for(auto res : output){ + it->data = res; + it++; + } for (z=1; z<=nz; ++z) { newPlane(); - i = 0 ; +// i = 0 ; cz = gridmin.z() + mstep * z ; - for (int y = 0 ; y < ny ; ++y) + //positions.clear(); +// for (int y = 0 ; y < ny ; ++y) +// { +// cy = gridmin.y() + mstep * y ; +// for (int x = 0 ; x < nx ; ++x, ++i) +// { +// cx = gridmin.x() + mstep * x ; +// //positions[i].set(cx, cy, cz); +// } +// } + + positions[0].z() = cz; + getFieldValueAt(positions, output) ; + + // Copy back the data into planes. + auto it = P1; + for(auto res : output) { - cy = gridmin.y() + mstep * y ; - for (int x = 0 ; x < nx ; ++x, ++i) - { - cx = gridmin.x() + mstep * x ; - - Vec3d pos { cx, cy, cz } ; - double res = getFieldValueAt(pos) ; - (P1+i)->data = res ; - } + it->data = res; + it++; } unsigned int i=0; diff --git a/applications/plugins/SofaImplicitField/MarchingCube.h b/applications/plugins/SofaImplicitField/MarchingCube.h index 85e2e6d49f6..6872a220f0f 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.h +++ b/applications/plugins/SofaImplicitField/MarchingCube.h @@ -40,7 +40,7 @@ class MarchingCube public: void generateSurfaceMesh(const double isoval, const double mstep, const double invStep, const Vec3d& gridmin, const Vec3d& gridmax, - std::function field, + std::function &, std::vector &)> field, SeqCoord& tmpPoints, SeqTriangles& tmpTriangles); private: diff --git a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp index 311e7816dd3..89a049f4b75 100644 --- a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp @@ -107,7 +107,13 @@ void FieldToSurfaceMesh::updateMeshIfNeeded() tmpTriangles.clear(); marchingCube.generateSurfaceMesh(isoval, mstep, invStep, gridmin, gridmax, - [field](Vec3d& pos){return field->getValue(pos);}, + [field](std::vector& positions, std::vector& res){ + res.reserve(positions.size()); + for(auto& position : positions) + { + res.emplace_back(field->getValue(position)); + } + }, tmpPoints, tmpTriangles); /// Copy the surface to Sofa topology diff --git a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl index f6467f631ef..267b73e3fc3 100644 --- a/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl +++ b/applications/plugins/SofaImplicitField/components/mapping/ImplicitSurfaceMapping.inl @@ -119,21 +119,30 @@ void ImplicitSurfaceMapping::apply(const core::MechanicalParams * /*mpar type::BoundingBox bigBox {mGridMin.getValue(), mGridMax.getValue()}; box.intersection(bigBox); - auto fieldFunction = [&sortParticles, &r, &r2, &invStep](Vec3d& pos) -> double { - int index = helper::rfloor(pos.z()*invStep); + auto fieldFunction = [&sortParticles, &r, &r2, &invStep]( + std::vector& pos, std::vector& res) -> void { + + auto z = pos[0].z(); + int index = helper::rfloor(z*invStep); auto particlesIt = sortParticles.find(index); if(particlesIt==sortParticles.end()) - return 0.0; - - double sumd = 0.0; - for(auto& particle : (particlesIt->second)){ - double d2 = (pos - particle).norm2(); - if(d2 < r2){ - d2 /= r2; - sumd += (1 + (-4*d2*d2*d2 + 17*d2*d2 - 22*d2)/9); + return; + + int i = 0; + for(auto& position : pos ) + { + 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); + } } + res[i++] = sumd; } - return sumd; + return; }; auto triangles = helper::getWriteOnlyAccessor(d_seqTriangles); @@ -145,7 +154,6 @@ void ImplicitSurfaceMapping::apply(const core::MechanicalParams * /*mpar invStep, box.minBBox(), box.maxBBox(), fieldFunction, points.wref(), triangles.wref()); - } template From c4058a006e1ca2d1cfe8b7cc3c54fbd67225f0ff Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Sat, 4 Oct 2025 21:07:21 +0200 Subject: [PATCH 11/22] Remove hacks --- .../SofaImplicitField/MarchingCube.cpp | 22 ++++++++----------- .../initSofaImplicitField.cpp | 1 - 2 files changed, 9 insertions(+), 14 deletions(-) diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp index 283ab21d3dc..6ef438a755f 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.cpp +++ b/applications/plugins/SofaImplicitField/MarchingCube.cpp @@ -100,20 +100,16 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, { newPlane(); -// i = 0 ; cz = gridmin.z() + mstep * z ; - //positions.clear(); -// for (int y = 0 ; y < ny ; ++y) -// { -// cy = gridmin.y() + mstep * y ; -// for (int x = 0 ; x < nx ; ++x, ++i) -// { -// cx = gridmin.x() + mstep * x ; -// //positions[i].set(cx, cy, cz); -// } -// } - - positions[0].z() = cz; + for (int y = 0 ; y < ny ; ++y) + { + cy = gridmin.y() + mstep * y ; + for (int x = 0 ; x < nx ; ++x, ++i) + { + cx = gridmin.x() + mstep * x ; + positions[i].set(cx, cy, cz); + } + } getFieldValueAt(positions, output) ; // Copy back the data into planes. diff --git a/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp b/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp index 13854029885..9e3a2d51e6f 100644 --- a/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp +++ b/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp @@ -103,7 +103,6 @@ const char* getModuleDescription() void registerObjects(sofa::core::ObjectFactory* factory) { - std::cout << "===================================================================" << std::endl; sofa::component::geometry::_BottleField_::registerBottleField(factory); sofa::component::geometry::_sphericalfield_::registerSphericalField(factory); sofa::component::geometry::_StarShapedField_::registerStarShapedField(factory); From ffeb514cb4d33f14d2c2570d008a699cd89d7a09 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Tue, 7 Oct 2025 21:55:44 +0200 Subject: [PATCH 12/22] Restored version --- .../SofaImplicitField/MarchingCube.cpp | 104 ++++++++---------- .../plugins/SofaImplicitField/MarchingCube.h | 2 - 2 files changed, 44 insertions(+), 62 deletions(-) diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp index 6ef438a755f..4aef282aaf9 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.cpp +++ b/applications/plugins/SofaImplicitField/MarchingCube.cpp @@ -30,21 +30,6 @@ namespace sofaimplicitfield { -void MarchingCube::newPlane() -{ - CubeData c; - c.p[0] = -1; - c.p[1] = -1; - c.p[2] = -1; - c.data = 0; - typename std::vector::iterator P = P0; - P0 = P1; - P1 = P; - int n = planes.size()/2; - for (int i=0; i&, std::vector&)> getFieldValueAt, @@ -54,39 +39,44 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, int ny = floor((gridmax.y() - gridmin.y()) * invStep) + 1 ; int nz = floor((gridmax.z() - gridmin.z()) * invStep) + 1 ; - if( nz <= 0 || ny <= 0 || nx <= 0 ) + // Marching cubes only works for a grid size larger than two + if( nz < 2 || ny < 2 || nx < 2 ) return; double cx,cy,cz; - int x,y,z,i,mk; + int z,mk; const int *tri; - planes.resize(2*(nx)*(ny)); - P0 = planes.begin()+0; + // Creates two planes + CubeData c{{-1,-1,-1},0}; + planes.reserve(2*nx*ny); + for(size_t i=0;i positions; - std::vector output; - positions.resize(ny*nx); - output.resize(nx*ny); - for (int y = 0 ; y < ny ; ++y) + positions.resize(nx*ny); + cz = gridmin.z(); + for (int i=0, y = 0 ; y < ny ; ++y) { cy = gridmin.y() + mstep * y ; - for (int x = 0 ; x < nx ; ++x, ++i) + for (int x = 0 ; x < nx ; ++x) { cx = gridmin.x() + mstep * x ; - positions[i].set(cx, cy, cz ); + positions[i++].set(cx, cy, cz ); } } + + std::vector output; getFieldValueAt(positions, output) ; // Copy back the data into planes. @@ -98,18 +88,19 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, for (z=1; z<=nz; ++z) { - newPlane(); + std::swap(P0, P1); cz = gridmin.z() + mstep * z ; - for (int y = 0 ; y < ny ; ++y) + for (int i=0, y=0 ; y < ny ; ++y) { cy = gridmin.y() + mstep * y ; - for (int x = 0 ; x < nx ; ++x, ++i) + for (int x = 0 ; x < nx ; ++x) { cx = gridmin.x() + mstep * x ; - positions[i].set(cx, cy, cz); + positions[i++].set(cx, cy, cz); } } + output.clear(); getFieldValueAt(positions, output) ; // Copy back the data into planes. @@ -120,7 +111,6 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, it++; } - unsigned int 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.begin(); @@ -139,56 +129,50 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, 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)) + if (((P1+di)->data>isoval)^((P1+di-dx)->data>isoval)) { - (P1+i)->p[0] = addPoint(tmpPoints, 0, pos,gridmin, (P1+i)->data,(P1+i-dx)->data, mstep, isoval); + (P1+di)->p[0] = addPoint(tmpPoints, 0, pos,gridmin, (P1+di)->data,(P1+di-dx)->data, mstep, isoval); } - if (((P1+i)->data>isoval)^((P1+i-dy)->data>isoval)) + if (((P1+di)->data>isoval)^((P1+di-dy)->data>isoval)) { - (P1+i)->p[1] = addPoint(tmpPoints, 1, pos,gridmin,(P1+i)->data,(P1+i-dy)->data, mstep, isoval); + (P1+di)->p[1] = addPoint(tmpPoints, 1, pos,gridmin,(P1+di)->data,(P1+di-dy)->data, mstep, isoval); } - if (((P1+i)->data>isoval)^((P0+i)->data>isoval)) + if (((P1+di)->data>isoval)^((P0+di)->data>isoval)) { - (P1+i)->p[2] = addPoint(tmpPoints, 2, pos,gridmin,(P1+i)->data,(P0+i)->data, mstep, 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+i-dx-dy)->data > isoval) mk = 1; + if ((P0+di-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; + 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+i; + 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; } - ++i; + ++di; } } } diff --git a/applications/plugins/SofaImplicitField/MarchingCube.h b/applications/plugins/SofaImplicitField/MarchingCube.h index 6872a220f0f..394a6163b1c 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.h +++ b/applications/plugins/SofaImplicitField/MarchingCube.h @@ -44,8 +44,6 @@ class MarchingCube SeqCoord& tmpPoints, SeqTriangles& tmpTriangles); private: - void newPlane(); - /// For each cube, store the vertex indices on each 3 first edges, and the data value struct CubeData { From 20b2355b5837ab45db48a6d1a355522cf162ff55 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Tue, 7 Oct 2025 22:42:45 +0200 Subject: [PATCH 13/22] [SofaImplicitField] FIX Binding_ScalarField invalid "name" set When the name is passed by key,value it is not set. --- .../plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp index fc2138831e4..927a801b7a3 100644 --- a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp +++ b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp @@ -112,6 +112,7 @@ void moduleAddScalarField(py::module &m) { "positional argument='" + py::cast(args[0]) + "'."); } + ff->setName(py::cast(value)); } } return ff; From 13b8b44cebf20a45cdc4ca2f44a58a81727db46a Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Tue, 7 Oct 2025 23:54:53 +0200 Subject: [PATCH 14/22] [SofaImplicitField] MarchingCube factorization and cleaning. --- applications/plugins/SofaImplicitField/MarchingCube.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp index 4aef282aaf9..02af01cc853 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.cpp +++ b/applications/plugins/SofaImplicitField/MarchingCube.cpp @@ -64,7 +64,9 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, z = 0; std::vector positions; + std::vector output; positions.resize(nx*ny); + output.resize(nx*ny); cz = gridmin.z(); for (int i=0, y = 0 ; y < ny ; ++y) { @@ -75,8 +77,6 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, positions[i++].set(cx, cy, cz ); } } - - std::vector output; getFieldValueAt(positions, output) ; // Copy back the data into planes. @@ -100,7 +100,6 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, positions[i++].set(cx, cy, cz); } } - output.clear(); getFieldValueAt(positions, output) ; // Copy back the data into planes. From 34f215cc364a8240be2f4062ebeaf6122ec86f5f Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Tue, 7 Oct 2025 23:55:15 +0200 Subject: [PATCH 15/22] [SofaImplicitField] Clean FieldToSurfaceMesh --- .../components/engine/FieldToSurfaceMesh.cpp | 32 +++++++++---------- .../components/engine/FieldToSurfaceMesh.h | 10 +++--- 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp index 89a049f4b75..a6f4a9cc290 100644 --- a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp @@ -39,17 +39,19 @@ FieldToSurfaceMesh::FieldToSurfaceMesh() , 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_dridMax(initData(&d_dridMax, Vec3d(1,1,1),"max","Grid Max")) + , 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_dridMax}, [this](const sofa::core::DataTracker&) + addUpdateCallback("updateMesh", {&d_step, &d_IsoValue, &d_gridMin, &d_gridMax}, [this](const sofa::core::DataTracker&) { checkInputs(); - hasChanged=true; + updateMeshIfNeeded(); return core::objectmodel::ComponentState::Valid; - }, {}); + }, {&d_outPoints, &d_outTriangles}); + d_outPoints.setGroup("Output"); + d_outTriangles.setGroup("Output"); } FieldToSurfaceMesh::~FieldToSurfaceMesh() @@ -64,14 +66,17 @@ void FieldToSurfaceMesh::init() 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_dridMax.getValue()-d_gridMin.getValue() ; + auto length = d_gridMax.getValue()-d_gridMin.getValue() ; auto step = d_step.getValue(); // clamp the mStep value to avoid too large grids @@ -84,9 +89,6 @@ void FieldToSurfaceMesh::checkInputs(){ void FieldToSurfaceMesh::updateMeshIfNeeded() { - if(!hasChanged) - return; - sofa::helper::getWriteOnlyAccessor(d_outPoints).clear(); sofa::helper::getWriteOnlyAccessor(d_outTriangles).clear(); @@ -95,7 +97,7 @@ void FieldToSurfaceMesh::updateMeshIfNeeded() double invStep = 1.0/d_step.getValue(); Vec3d gridmin = d_gridMin.getValue() ; - Vec3d gridmax = d_dridMax.getValue() ; + Vec3d gridmax = d_gridMax.getValue() ; auto field = l_field.get(); @@ -108,10 +110,10 @@ void FieldToSurfaceMesh::updateMeshIfNeeded() marchingCube.generateSurfaceMesh(isoval, mstep, invStep, gridmin, gridmax, [field](std::vector& positions, std::vector& res){ - res.reserve(positions.size()); + int i=0; for(auto& position : positions) { - res.emplace_back(field->getValue(position)); + res[i++]=field->getValue(position); } }, tmpPoints, tmpTriangles); @@ -135,12 +137,8 @@ void FieldToSurfaceMesh::draw(const VisualParams* vparams) if(!d_debugDraw.getValue()) return; - updateMeshIfNeeded(); - auto drawTool = vparams->drawTool(); - drawTool->drawBoundingBox(d_gridMin.getValue(), d_dridMax.getValue()) ; - sofa::helper::ReadAccessor< Data > x = d_outPoints; sofa::helper::ReadAccessor< Data > triangles = d_outTriangles; drawTool->setLightingEnabled(true); diff --git a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h index 17916561cca..c4dbc73882e 100644 --- a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h @@ -58,9 +58,9 @@ class FieldToSurfaceMesh : public BaseObject 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_dridMax.getValue(); } - void setGridMax(const Vec3d& val) { d_dridMax.setValue(val); } - void setGridMax(double x, double y, double z) { d_dridMax.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 d_IsoValue; Data< Vec3d > d_gridMin; - Data< Vec3d > d_dridMax; + Data< Vec3d > d_gridMax; /// Output Data d_outPoints; @@ -82,6 +82,8 @@ class FieldToSurfaceMesh : public BaseObject virtual ~FieldToSurfaceMesh() ; private: + void computeBBox(const core::ExecParams* /* params */, bool /*onlyVisible*/=false) override; + void checkInputs(); void generateSurfaceMesh(double isoval, double mstep, double invStep, From e80bd2bf56770aa48c18e4a27ed06deb03774d39 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Tue, 7 Oct 2025 23:58:01 +0200 Subject: [PATCH 16/22] [SofaImplicitField] Add few field function and operator in the xshape python package --- .../examples/python/primitives.py | 14 -------- .../examples/python/xshape/__init__.py | 0 .../examples/python/xshape/operators.py | 35 +++++++++++++++++++ .../examples/python/xshape/primitives.py | 34 ++++++++++++++++++ .../examples/python/xshape/transforms.py | 15 ++++++++ 5 files changed, 84 insertions(+), 14 deletions(-) delete mode 100644 applications/plugins/SofaImplicitField/examples/python/primitives.py create mode 100644 applications/plugins/SofaImplicitField/examples/python/xshape/__init__.py create mode 100644 applications/plugins/SofaImplicitField/examples/python/xshape/operators.py create mode 100644 applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py create mode 100644 applications/plugins/SofaImplicitField/examples/python/xshape/transforms.py diff --git a/applications/plugins/SofaImplicitField/examples/python/primitives.py b/applications/plugins/SofaImplicitField/examples/python/primitives.py deleted file mode 100644 index 174e9ddbd7d..00000000000 --- a/applications/plugins/SofaImplicitField/examples/python/primitives.py +++ /dev/null @@ -1,14 +0,0 @@ - -from SofaImplicitField import ScalarField -import numpy - -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.sqrt( numpy.sum((self.center.value - numpy.array([x,y,z]))**2) ) - self.radius.value \ No newline at end of file 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..87f8a34fff2 --- /dev/null +++ b/applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py @@ -0,0 +1,34 @@ +"""Distance field function + +Sources: + https://iquilezles.org/articles/distfunctions/ +""" +from SofaImplicitField import ScalarField +import numpy + +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 + +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 From 044a11db353247e7f1d2f9b96c66017bd39679eb Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Tue, 7 Oct 2025 23:58:30 +0200 Subject: [PATCH 17/22] [SofaImplicitField] Add an example of use of the FieldToSurfaceMesh component --- .../example-mesh-extraction-from-implicit.py | 31 ++++++++++++++++--- 1 file changed, 27 insertions(+), 4 deletions(-) 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 index 895a7191165..02d506e2066 100644 --- a/applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py +++ b/applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py @@ -1,6 +1,18 @@ import Sofa -from primitives import Sphere +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++ @@ -8,11 +20,22 @@ def createScene(root : Sofa.Core.Node): """ root.addObject("RequiredPlugin", name="SofaImplicitField") + root.addObject(DrawController()) + ########################### Fields ################## root.addChild("Fields") - f1 = root.Fields.addObject(Sphere(name="field1", center=[0,0,0])) - f2 = root.Fields.addObject("SphericalField", name="field2", center=[2,0,0]) + 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", @@ -21,7 +44,7 @@ def createScene(root : Sofa.Core.Node): m2 = root.Meshing.addObject("FieldToSurfaceMesh", name="polygonizer2", field=f2.linkpath, min=[1,-1,-1], max=[3,1,1], - step=0.01) + step=0.07) ########################### Fields ################## root.addChild("Visual") From a870250a66c74816f9807329853d888bf474cbb7 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 8 Oct 2025 11:01:41 +0200 Subject: [PATCH 18/22] [SofaImplicitField] Fix memory allocation bug. --- .../plugins/SofaImplicitField/MarchingCube.cpp | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp index 02af01cc853..071407b2376 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.cpp +++ b/applications/plugins/SofaImplicitField/MarchingCube.cpp @@ -49,12 +49,12 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, // Creates two planes CubeData c{{-1,-1,-1},0}; - planes.reserve(2*nx*ny); + planes.resize(2*nx*ny); for(size_t i=0;idata = res; it++; @@ -88,8 +88,6 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, for (z=1; z<=nz; ++z) { - std::swap(P0, P1); - cz = gridmin.z() + mstep * z ; for (int i=0, y=0 ; y < ny ; ++y) { @@ -174,6 +172,7 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, ++di; } } + std::swap(P0, P1); } } From d664bd9573ce0cd8c87070da59dd91f843094a68 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 8 Oct 2025 11:02:13 +0200 Subject: [PATCH 19/22] [SofaImplicitField] FIX the rendering of normal that was inverted. --- .../components/engine/FieldToSurfaceMesh.cpp | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp index a6f4a9cc290..e42e18467a9 100644 --- a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp @@ -153,11 +153,7 @@ void FieldToSurfaceMesh::draw(const VisualParams* vparams) Vec3d pb = (0.9*x[b]+0.1*center) ; Vec3d pc = (0.9*x[c]+0.1*center) ; - Vec3d a1 = x[c]-x[b] ; - Vec3d a2 = x[a]-x[b] ; - - vparams->drawTool()->drawTriangles({pa,pb,pc}, - a1.cross(a2), + vparams->drawTool()->drawTriangles({pb,pa,pc}, type::RGBAColor(0.0,0.0,1.0,1.0)); } From 66cfe3fb8cb477d7d582e2e4ac02d2f209da3722 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 8 Oct 2025 11:15:21 +0200 Subject: [PATCH 20/22] [SofaImplicitField] Factorize the filling code in a lambda. --- .../SofaImplicitField/MarchingCube.cpp | 56 +++++++------------ 1 file changed, 21 insertions(+), 35 deletions(-) diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp index 071407b2376..4a7d0098b20 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.cpp +++ b/applications/plugins/SofaImplicitField/MarchingCube.cpp @@ -63,50 +63,36 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, z = 0; - std::vector positions; - std::vector output; - positions.resize(nx*ny); - output.resize(nx*ny); - cz = gridmin.z(); - for (int i=0, y = 0 ; y < ny ; ++y) + auto fillPlane = [getFieldValueAt](std::vector &positions, std::vector& output, + double mstep, double gridmin_y, double gridmin_x, int ny, int nx, float cz, + std::vector::iterator itDestPlane) { - cy = gridmin.y() + mstep * y ; - for (int x = 0 ; x < nx ; ++x) + for (int i=0, y = 0 ; y < ny ; ++y) { - cx = gridmin.x() + mstep * x ; - positions[i++].set(cx, cy, cz ); - } - } - getFieldValueAt(positions, output) ; - - // Copy back the data into planes. - auto it = P0; - for(auto res : output){ - it->data = res; - it++; - } - - for (z=1; z<=nz; ++z) - { - cz = gridmin.z() + mstep * z ; - for (int i=0, y=0 ; y < ny ; ++y) - { - cy = gridmin.y() + mstep * y ; + double cy = gridmin_y + mstep * y ; for (int x = 0 ; x < nx ; ++x) { - cx = gridmin.x() + mstep * x ; - positions[i++].set(cx, cy, cz); + double cx = gridmin_x + mstep * x ; + positions[i++].set(cx, cy, cz ); } } getFieldValueAt(positions, output) ; - // Copy back the data into planes. - auto it = P1; - for(auto res : output) - { - it->data = res; - it++; + 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}; From e966ddd930ca4db59d04d1dd318a245d7c2d93c5 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 8 Oct 2025 21:20:50 +0200 Subject: [PATCH 21/22] Add vectorized version --- .../SofaImplicitField/MarchingCube.cpp | 1 - .../components/engine/FieldToSurfaceMesh.cpp | 9 +-- .../components/geometry/ScalarField.cpp | 11 +++ .../components/geometry/ScalarField.h | 3 + .../examples/python/xshape/primitives.py | 9 +++ .../python/src/Binding_ScalarField.cpp | 70 +++++++++++++++++++ 6 files changed, 96 insertions(+), 7 deletions(-) diff --git a/applications/plugins/SofaImplicitField/MarchingCube.cpp b/applications/plugins/SofaImplicitField/MarchingCube.cpp index 4a7d0098b20..a9b74985655 100644 --- a/applications/plugins/SofaImplicitField/MarchingCube.cpp +++ b/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; diff --git a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp index e42e18467a9..fbf12bdbf84 100644 --- a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp +++ b/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& positions, std::vector& res){ - int i=0; - for(auto& position : positions) - { - res[i++]=field->getValue(position); - } + [field](const std::vector& positions, std::vector& res){ + field->getValues(positions, res); }, tmpPoints, tmpTriangles); diff --git a/applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp b/applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp index 9c0323bbfd8..2eae46b226b 100644 --- a/applications/plugins/SofaImplicitField/components/geometry/ScalarField.cpp +++ b/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& 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 bebd099e6f0..aa27b5dca04 100644 --- a/applications/plugins/SofaImplicitField/components/geometry/ScalarField.h +++ b/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& 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/examples/python/xshape/primitives.py b/applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py index 87f8a34fff2..fe203d79654 100644 --- a/applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py +++ b/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) diff --git a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp index 927a801b7a3..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); From 4f4e66dc095dfa431fc7087fe0c84f9881c076c6 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Thu, 7 May 2026 16:41:41 +0200 Subject: [PATCH 22/22] FIXs --- .../examples/python/example-mesh-extraction-from-implicit.py | 2 +- .../SofaImplicitField/examples/python/xshape/primitives.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) 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 index 7ad2d606c5f..4f3c9d07610 100644 --- a/applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py +++ b/applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py @@ -18,7 +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 ################## diff --git a/applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py b/applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py index fe203d79654..1321ff3ad5b 100644 --- a/applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py +++ b/applications/plugins/SofaImplicitField/examples/python/xshape/primitives.py @@ -5,7 +5,6 @@ """ from SofaImplicitField import ScalarField import numpy -import numba class Sphere(ScalarField): def __init__(self, *args, **kwargs):