From 51e310f079d44886b4d9ceb0e90311922a358bd3 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 25 Jun 2025 00:02:37 +0200 Subject: [PATCH 1/8] [SofaImplicitField] Add python binding to SofaImplicitField --- .../plugins/SofaImplicitField/CMakeLists.txt | 4 + .../examples/python/python-scalarfield.py | 15 ++++ .../SofaImplicitField/python/CMakeLists.txt | 26 ++++++ .../python/src/Binding_ScalarField.cpp | 80 +++++++++++++++++++ .../python/src/Binding_ScalarField.h | 29 +++++++ .../python/src/Module_SofaImplicitField.cpp | 38 +++++++++ 6 files changed, 192 insertions(+) create mode 100644 applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py create mode 100644 applications/plugins/SofaImplicitField/python/CMakeLists.txt create mode 100644 applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp create mode 100644 applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.h create mode 100644 applications/plugins/SofaImplicitField/python/src/Module_SofaImplicitField.cpp diff --git a/applications/plugins/SofaImplicitField/CMakeLists.txt b/applications/plugins/SofaImplicitField/CMakeLists.txt index f7cfab05c64..35d08e7f7b1 100644 --- a/applications/plugins/SofaImplicitField/CMakeLists.txt +++ b/applications/plugins/SofaImplicitField/CMakeLists.txt @@ -47,6 +47,10 @@ endif() add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES} ${EXTRA_FILES}) target_link_libraries(${PROJECT_NAME} PRIVATE Sofa.Component.Topology.Container.Constant) +find_package(SofaPython3 REQUIRED) +if (SofaPython3_FOUND) + add_subdirectory(python) +endif() ## Install rules for the library and headers; CMake package configurations files sofa_create_package_with_targets( diff --git a/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py b/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py new file mode 100644 index 00000000000..ba044b8abda --- /dev/null +++ b/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py @@ -0,0 +1,15 @@ +import Sofa +from SofaImplicitField import ScalarField + +class Sphere(ScalarField): + def __init__(self, *args, **kwargs): + ScalarField.__init__(self, *args, **kwargs) + self.addNewData("center") + + + def getValue(self, x, y, z): + return 123.456 + +def createScene(root): + root.addObject(MyField("field")) + print( "Field value is: ", root.field.getValue(0.0,0.0,0.0) ) \ No newline at end of file diff --git a/applications/plugins/SofaImplicitField/python/CMakeLists.txt b/applications/plugins/SofaImplicitField/python/CMakeLists.txt new file mode 100644 index 00000000000..bd66bd9546b --- /dev/null +++ b/applications/plugins/SofaImplicitField/python/CMakeLists.txt @@ -0,0 +1,26 @@ +project(SofaImplicitField.Python) + +set(SOURCE_FILES + ${CMAKE_CURRENT_SOURCE_DIR}/src/Binding_ScalarField.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/Module_SofaImplicitField.cpp +) + +set(HEADER_FILES + ${CMAKE_CURRENT_SOURCE_DIR}/src/Binding_ScalarField.h +) + +if (NOT TARGET SofaPython3::Plugin) + find_package(SofaPython3 REQUIRED COMPONENTS Plugin Bindings.Sofa.Core) +endif() + +sofa_find_package(SofaImplicitField REQUIRED) + +SP3_add_python_module( + TARGET ${PROJECT_NAME} + PACKAGE SofaImplicitField.Python + MODULE SofaImplicitField + DESTINATION . + SOURCES ${SOURCE_FILES} + HEADERS ${HEADER_FILES} + DEPENDS SofaImplicitField SofaPython3::Plugin SofaPython3::Bindings.Sofa.Core +) diff --git a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp new file mode 100644 index 00000000000..21df8bb41d6 --- /dev/null +++ b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp @@ -0,0 +1,80 @@ +#include + +#include +#include +#include +#include + +#include "Binding_ScalarField.h" +/****************************************************************************** +* SofaImplicitField plugin * +* (c) 2024 CNRS, University of Lille, INRIA * +* * +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +/// Makes an alias for the pybind11 namespace to increase readability. +namespace py { using namespace pybind11; } + +namespace sofaimplicitfield { +using namespace sofapython3; +using sofa::component::geometry::ScalarField; +using sofa::core::objectmodel::BaseObject; +using sofa::type::Vec3d; + +class ScalarField_Trampoline : public ScalarField { +public: + SOFA_CLASS(ScalarField_Trampoline, ScalarField); + + double getValue(Vec3d& pos, int& domain) override{ + SOFA_UNUSED(domain); + PythonEnvironment::gil acquire; + + PYBIND11_OVERLOAD_PURE(double, ScalarField, getValue, pos.x(), pos.y(), pos.z()); + } +}; + +void moduleAddScalarField(py::module &m) { + py::class_> f(m, "ScalarField", py::dynamic_attr(), ""); + + f.def(py::init([](py::args &args, py::kwargs &kwargs) { + auto ff = sofa::core::sptr (new ScalarField_Trampoline()); + + ff->f_listening.setValue(true); + + if (args.size() == 1) ff->setName(py::cast(args[0])); + + py::object cc = py::cast(ff); + for (auto kv : kwargs) { + std::string key = py::cast(kv.first); + py::object value = py::reinterpret_borrow(kv.second); + if (key == "name") { + if (args.size() != 0) { + throw py::type_error("The name is set twice as a " + "named argument='" + py::cast(value) + "' and as a" + "positional argument='" + + py::cast(args[0]) + "'."); + } + } + } + return ff; + })); + + m.def("getValue", &ScalarField_Trampoline::getValue); +} + +} diff --git a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.h b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.h new file mode 100644 index 00000000000..8ab31e32945 --- /dev/null +++ b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.h @@ -0,0 +1,29 @@ +/****************************************************************************** +* SofaImplicitField plugin * +* (c) 2021 CNRS, University of Lille, INRIA * +* * +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#pragma once + +#include + +namespace sofaimplicitfield { + +void moduleAddScalarField(pybind11::module &m); + +} diff --git a/applications/plugins/SofaImplicitField/python/src/Module_SofaImplicitField.cpp b/applications/plugins/SofaImplicitField/python/src/Module_SofaImplicitField.cpp new file mode 100644 index 00000000000..bd7e9480856 --- /dev/null +++ b/applications/plugins/SofaImplicitField/python/src/Module_SofaImplicitField.cpp @@ -0,0 +1,38 @@ +/****************************************************************************** +* SofaImplicitField plugin * +* (c) 2024 CNRS, University of Lille, INRIA * +* * +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#include +namespace py = pybind11; + +#include "Binding_ScalarField.h" + +namespace sofaimplicitfield +{ + +PYBIND11_MODULE(SofaImplicitField, m) { + m.doc() = R"doc( + Implement scalar field representation in python + )doc"; + + moduleAddScalarField(m); +} + +} + From 34dd7f301a166836b5ad6f86f37afc86d61e95a7 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 25 Jun 2025 00:17:53 +0200 Subject: [PATCH 2/8] Update example. --- .../examples/python/python-scalarfield.py | 23 +++++++++++++++---- .../python/src/Binding_ScalarField.cpp | 16 ++++++------- 2 files changed, 26 insertions(+), 13 deletions(-) diff --git a/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py b/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py index ba044b8abda..d97e36a5814 100644 --- a/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py +++ b/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py @@ -1,15 +1,28 @@ import Sofa from SofaImplicitField import ScalarField +import numpy class Sphere(ScalarField): def __init__(self, *args, **kwargs): ScalarField.__init__(self, *args, **kwargs) - self.addNewData("center") - + + 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 123.456 + 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(MyField("field")) - print( "Field value is: ", root.field.getValue(0.0,0.0,0.0) ) \ No newline at end of file + root.addObject(Sphere("field")) + root.addObject(FieldController(target=root.field)) diff --git a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp index 21df8bb41d6..9e2243bd5ab 100644 --- a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp +++ b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp @@ -1,11 +1,3 @@ -#include - -#include -#include -#include -#include - -#include "Binding_ScalarField.h" /****************************************************************************** * SofaImplicitField plugin * * (c) 2024 CNRS, University of Lille, INRIA * @@ -25,6 +17,14 @@ ******************************************************************************* * Contact information: contact@sofa-framework.org * ******************************************************************************/ +#include + +#include +#include +#include +#include + +#include "Binding_ScalarField.h" /// Makes an alias for the pybind11 namespace to increase readability. namespace py { using namespace pybind11; } From f815374b21e93cd31b8484b2cc55f38939dcc36e Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Thu, 26 Jun 2025 21:49:25 +0200 Subject: [PATCH 3/8] [Core] Add ObjectFactory in sofa/core/fwd --- Sofa/framework/Core/src/sofa/core/fwd.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sofa/framework/Core/src/sofa/core/fwd.h b/Sofa/framework/Core/src/sofa/core/fwd.h index 622a30e1394..dc2c5109e47 100644 --- a/Sofa/framework/Core/src/sofa/core/fwd.h +++ b/Sofa/framework/Core/src/sofa/core/fwd.h @@ -30,7 +30,7 @@ namespace sofa::helper::visual { class DrawTool; } namespace sofa::core { - +class ObjectFactory; class BaseState; class BaseMapping; class BehaviorModel; From f9c69a3bfb95ffcef1b19e586518c6b3d5861065 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Thu, 26 Jun 2025 21:50:47 +0200 Subject: [PATCH 4/8] [Core] Add a getRevisionCounter() to hide details of who stores this counter. This is then used to connected and track component revisions. --- Sofa/framework/Core/src/sofa/core/objectmodel/Base.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sofa/framework/Core/src/sofa/core/objectmodel/Base.h b/Sofa/framework/Core/src/sofa/core/objectmodel/Base.h index ae3b2285b6c..330a3d7cc4b 100644 --- a/Sofa/framework/Core/src/sofa/core/objectmodel/Base.h +++ b/Sofa/framework/Core/src/sofa/core/objectmodel/Base.h @@ -352,7 +352,7 @@ class SOFA_CORE_API Base : public IntrusiveObject /// @name componentstate /// Methods related to component state /// @{ - + int getRevisionCounter() const { return d_componentState.getCounter(); } ComponentState getComponentState() const { return d_componentState.getValue() ; } bool isComponentStateValid() const { return d_componentState.getValue() == ComponentState::Valid; } bool isComponentStateInvalid() const { return d_componentState.getValue() == ComponentState::Invalid; } From e4a5250972ae66f3d68c49873bb8bf843c18afdc Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Thu, 26 Jun 2025 21:51:43 +0200 Subject: [PATCH 5/8] [Component.Visual] Add a new update callback in VisualModelImpl Currently, when the triangles set is changed, the geometry is not modified. --- .../Visual/src/sofa/component/visual/VisualModelImpl.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/Sofa/Component/Visual/src/sofa/component/visual/VisualModelImpl.cpp b/Sofa/Component/Visual/src/sofa/component/visual/VisualModelImpl.cpp index fc93fac1094..1b2ce730267 100644 --- a/Sofa/Component/Visual/src/sofa/component/visual/VisualModelImpl.cpp +++ b/Sofa/Component/Visual/src/sofa/component/visual/VisualModelImpl.cpp @@ -166,10 +166,15 @@ VisualModelImpl::VisualModelImpl() //const std::string &name, std::string filena // add one identity matrix xforms.resize(1); + addUpdateCallback("updateGeometry", {&d_triangles}, + [&](const core::DataTracker&) -> sofa::core::objectmodel::ComponentState { + modified=true; + return sofa::core::objectmodel::ComponentState::Loading; + }, {&d_componentState}); + addUpdateCallback("updateTextures", { &d_texturename }, - [&](const core::DataTracker& tracker) -> sofa::core::objectmodel::ComponentState + [&](const core::DataTracker&) -> sofa::core::objectmodel::ComponentState { - SOFA_UNUSED(tracker); m_textureChanged = true; return sofa::core::objectmodel::ComponentState::Loading; }, { &d_componentState }); From 369c817107bc1522ae80b655738da9e2bab4ecb2 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Thu, 26 Jun 2025 21:55:14 +0200 Subject: [PATCH 6/8] =?UTF-8?q?[SofaImplicitField]=C2=A0Add=20component=20?= =?UTF-8?q?FieldToSurfaceMesh=20to=20produce=20a=20triangulation=20of=20a?= =?UTF-8?q?=20field.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The computation can be made asynchronous, so multiple rasterisation can be done in parallel. Depending on the connected component this may not work because the connected component may not handle dynamic geometry. --- .../plugins/SofaImplicitField/CMakeLists.txt | 5 + .../components/engine/FieldToSurfaceMesh.cpp | 356 ++++++++++++++++++ .../components/engine/FieldToSurfaceMesh.h | 135 +++++++ .../components/engine/RayMarching.cpp | 0 .../components/engine/RayMarching.h | 0 .../components/geometry/ScalarField.cpp | 5 + applications/plugins/SofaImplicitField/fwd.h | 33 ++ .../initSofaImplicitField.cpp | 4 + 8 files changed, 538 insertions(+) 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/components/engine/RayMarching.cpp create mode 100644 applications/plugins/SofaImplicitField/components/engine/RayMarching.h create mode 100644 applications/plugins/SofaImplicitField/fwd.h diff --git a/applications/plugins/SofaImplicitField/CMakeLists.txt b/applications/plugins/SofaImplicitField/CMakeLists.txt index 35d08e7f7b1..a3d74a87fa4 100644 --- a/applications/plugins/SofaImplicitField/CMakeLists.txt +++ b/applications/plugins/SofaImplicitField/CMakeLists.txt @@ -5,6 +5,7 @@ sofa_find_package(Sofa.Component.Topology.Container.Constant REQUIRED) set(HEADER_FILES config.h.in + fwd.h initSofaImplicitField.h # This is backward compatibility @@ -12,6 +13,8 @@ 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/engine/RayMarching.h components/geometry/BottleField.h components/geometry/DiscreteGridField.h components/geometry/SphericalField.h @@ -28,6 +31,8 @@ set(SOURCE_FILES deprecated/SphereSurface.cpp deprecated/InterpolatedImplicitSurface.cpp + components/engine/FieldToSurfaceMesh.cpp + components/engine/RayMarching.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..a7a4ac9cf5b --- /dev/null +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.cpp @@ -0,0 +1,356 @@ +/****************************************************************************** +* 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" + +#include + +namespace sofa::component::geometry +{ + +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_out_points(initData(&d_out_points, "outputPoints", "position of the tiangles vertex")) + , d_out_triangles(initData(&d_out_triangles, "outputTriangles", "list of triangles")) + , d_debugDraw(initData(&d_debugDraw,true, "debugDraw","Display the extracted surface")) + , d_doAsync(initData(&d_doAsync,false, "doAsync","Extract the surface in an asynchronous thread")) +{ + addUpdateCallback("updateMesh", {&mStep, &mIsoValue, &mGridMin, &mGridMax}, [this](const core::DataTracker&){ + checkInputs(); + + std::cout << " UPDATE MESH " << getRevisionCounter() << std::endl; + 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; + } + 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(){ + // Check if we are computing, if so, wait termination + if(workInProgress.load()) + return; + + // Check that the component state has not its revision counter changed. + if( !workFinished.load() && + (l_field->getRevisionCounter() != lastGenerationFieldCounter || getRevisionCounter() != lastGenerationCounter) ) + { + std::cout << "STARTING A MESH " << getName() << std::endl; + // First clean the two buffers. + sofa::helper::getWriteOnlyAccessor(d_out_points).clear(); + sofa::helper::getWriteOnlyAccessor(d_out_triangles).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(); + results = std::async(std::launch::async, [isoval, mstep, invStep, gridmin, gridmax, field, + this](){ + workInProgress.store(true); + generateSurfaceMesh(isoval, mstep, invStep, gridmin, gridmax, field); + workFinished.store(true); + workInProgress.store(false); + return true; + }); + + if(!d_doAsync.getValue()) + results.wait(); + } + + if(workFinished.load()){ + /// Copy the surface to Sofa topology + d_out_points.setValue(tmpPoints); + d_out_triangles.setValue(tmpTriangles); + + tmpPoints.clear(); + tmpTriangles.clear(); + + lastGenerationFieldCounter = l_field->getRevisionCounter(); + lastGenerationCounter = getRevisionCounter(); + workFinished.store(false); + //results.wait(); + //results.get(); + std::cout << "FINISHING A MESH " << getName() + << " " << d_out_points.getValue().size() << " " << d_out_triangles.getValue().size() << std::endl; + + return; + } +} + +void FieldToSurfaceMesh::draw(const VisualParams* vparams) +{ + if(isComponentStateInvalid()) + return; + + if(!d_debugDraw.getValue()) + return; + + updateMeshIfNeeded(); + + auto dt = vparams->drawTool(); + + dt->drawBoundingBox(mGridMin.getValue(), mGridMax.getValue()) ; + + sofa::helper::ReadAccessor< Data > x = d_out_points; + sofa::helper::ReadAccessor< Data > triangles = d_out_triangles; + dt->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()>10000){ + dt->drawPoints(x, 1.0, type::RGBAColor(1.0,1.0,0.0,0.2)) ; + }else{ + dt->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; + + std::cout << "1" << std::endl; + 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; + + + std::cout << "2" << std::endl; + 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; + if (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())<0) + { + /*std::cout << " mk=0x"<()); +} + +} diff --git a/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h new file mode 100644 index 00000000000..28d57a40efe --- /dev/null +++ b/applications/plugins/SofaImplicitField/components/engine/FieldToSurfaceMesh.h @@ -0,0 +1,135 @@ +/****************************************************************************** +* 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 sofa::component::geometry +{ + +typedef sofa::core::topology::BaseMeshTopology::SeqTriangles SeqTriangles; +typedef sofa::core::topology::BaseMeshTopology::Triangle Triangle; +typedef sofa::type::vector VecCoord; + +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_out_points; + Data d_out_triangles; + Data d_debugDraw; + Data d_doAsync; + + 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(); + void generateSurfaceMesh(double isoval, double mstep, double invStep, + Vec3d gridmin, Vec3d gridmax, + sofa::component::geometry::ScalarField*); + void updateMeshIfNeeded(); + +protected: + FieldToSurfaceMesh() ; + virtual ~FieldToSurfaceMesh() ; + + int lastGenerationCounter {-1}; + int lastGenerationFieldCounter {-1}; + +private: + void checkInputs(); + std::atomic workInProgress; + std::atomic workFinished; + std::future results; + VecCoord tmpPoints; + SeqTriangles tmpTriangles; +}; + +} + diff --git a/applications/plugins/SofaImplicitField/components/engine/RayMarching.cpp b/applications/plugins/SofaImplicitField/components/engine/RayMarching.cpp new file mode 100644 index 00000000000..e69de29bb2d diff --git a/applications/plugins/SofaImplicitField/components/engine/RayMarching.h b/applications/plugins/SofaImplicitField/components/engine/RayMarching.h new file mode 100644 index 00000000000..e69de29bb2d 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/fwd.h b/applications/plugins/SofaImplicitField/fwd.h new file mode 100644 index 00000000000..a026ff9ff3c --- /dev/null +++ b/applications/plugins/SofaImplicitField/fwd.h @@ -0,0 +1,33 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 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 + +namespace sofa::component::geometry +{ + class FieldToSurfaceMesh; + + template void registerToFactory(sofa::core::ObjectFactory* factory); +} + diff --git a/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp b/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp index ca5f0f5732d..f9ee30d23e2 100644 --- a/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp +++ b/applications/plugins/SofaImplicitField/initSofaImplicitField.cpp @@ -20,12 +20,15 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #include +#include #include #include #include using sofa::helper::system::PluginManager ; +using namespace sofa::component::geometry; + namespace sofa::component::geometry::_BottleField_ { extern void registerBottleField(sofa::core::ObjectFactory* factory); @@ -106,6 +109,7 @@ void registerObjects(sofa::core::ObjectFactory* factory) sofa::component::mapping::registerImplicitSurfaceMapping(factory); sofa::component::container::registerInterpolatedImplicitSurface(factory); sofa::component::geometry::_discretegrid_::registerDiscreteGridField(factory); + sofa::component::geometry::registerToFactory(factory); } } /// sofaimplicitfield From 05963cf34b25f6b9d1fa726ec501e14b66a8ca5c Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Thu, 26 Jun 2025 21:56:07 +0200 Subject: [PATCH 7/8] =?UTF-8?q?[SofaImplicitFIeld]=C2=A0Add=20examples=20o?= =?UTF-8?q?f=20use=20for=20component=20FieldToSurfaceMesh?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../example-mesh-extraction-from-implicit.py | 23 +++++++++++++++++ .../examples/python/python-scalarfield.py | 23 ++++++++--------- .../examples/python/shapes.py | 25 +++++++++++++++++++ 3 files changed, 59 insertions(+), 12 deletions(-) create mode 100644 applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py create mode 100644 applications/plugins/SofaImplicitField/examples/python/shapes.py 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..d297f84683a --- /dev/null +++ b/applications/plugins/SofaImplicitField/examples/python/example-mesh-extraction-from-implicit.py @@ -0,0 +1,23 @@ +import Sofa +from shapes import Sphere + +def createScene(root : Sofa.Core.Node): + root.addObject("RequiredPlugin", name="SofaImplicitField") + + s = root.addObject(Sphere(name="field1", center=[0,0,0])) + root.addObject("SphericalField", name="field2", center=[2,0,0]) + + root.addChild("Visual") + + root.Visual.addObject("FieldToSurfaceMesh", name="polygonizer1", + field=s.linkpath, min=[-1,-1,-1], max=[1,1,1], + isoValue="0.0", step="0.1",doAsync=True) + + root.Visual.addObject("FieldToSurfaceMesh", name="polygonizer2", + field=root.field2.linkpath, min=[1,-1,-1], max=[3,1,1], + isoValue="0.0", step="0.01",doAsync=True) + + root.Visual.addObject("OglModel", name="renderer", + position=root.Visual.polygonizer2.outputPoints.linkpath, + triangles=root.Visual.polygonizer2.outputTriangles.linkpath) + \ No newline at end of file diff --git a/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py b/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py index d97e36a5814..4de1988e932 100644 --- a/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py +++ b/applications/plugins/SofaImplicitField/examples/python/python-scalarfield.py @@ -1,16 +1,5 @@ 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 +from Shapes import Sphere class FieldController(Sofa.Core.Controller): def __init__(self, *args, **kwargs): @@ -23,6 +12,16 @@ def onAnimateEndEvent(self, event): 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 gradient at 1.0,0.0,0.0 is: ", self.field.getGradient([2.0,0.0,0.0]) ) + print("Field hessian at 1.0,0.0,0.0 is: ", self.field.getHessian([2.0,0.0,0.0]) ) + def createScene(root): root.addObject(Sphere("field")) root.addObject(FieldController(target=root.field)) + + root.addChild("Visual") + root.Visual.addObject("OglModel", name="renderer") + root.Visual.addObject("ImplicitSurfaceMapping", name="polygonizer", + input=root., output=root.Visual.renderer.linkpath, + isoValue="0.5", radius="0.75", step="0.25") + \ No newline at end of file diff --git a/applications/plugins/SofaImplicitField/examples/python/shapes.py b/applications/plugins/SofaImplicitField/examples/python/shapes.py new file mode 100644 index 00000000000..ce8593cee88 --- /dev/null +++ b/applications/plugins/SofaImplicitField/examples/python/shapes.py @@ -0,0 +1,25 @@ +import Sofa +from SofaImplicitField import ScalarField +import numpy +import jax +import jax.numpy as jnp + +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") + + self._center = jnp.array(self.center.value) + + self.fast_getValue = jax.jit(self.getValueFct) + + #def getValue(self, x, y, z): + # return numpy.sqrt( numpy.sum((self.center.value - numpy.array([x,y,z]))**2) ) - self.radius.value + + def getValueFct(self, x, y, z): + return jnp.sqrt( jnp.sum((self._center - jnp.array([x,y,z]))**2) ) - self.radius.value + + def getValue(self, x, y, z): + return self.fast_getValue(x, y, z) \ No newline at end of file From 6c8c226e7c20ac3926289c47ebf9b83acdecfe2d Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Thu, 26 Jun 2025 21:56:41 +0200 Subject: [PATCH 8/8] [SofaImplicitField.Binding] add the binding of the gradien and the hessian --- .../python/src/Binding_ScalarField.cpp | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp index 9e2243bd5ab..7a52fad43dc 100644 --- a/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp +++ b/applications/plugins/SofaImplicitField/python/src/Binding_ScalarField.cpp @@ -34,6 +34,7 @@ using namespace sofapython3; using sofa::component::geometry::ScalarField; using sofa::core::objectmodel::BaseObject; using sofa::type::Vec3d; +using sofa::type::Mat3x3; class ScalarField_Trampoline : public ScalarField { public: @@ -42,9 +43,21 @@ class ScalarField_Trampoline : public ScalarField { double getValue(Vec3d& pos, int& domain) override{ SOFA_UNUSED(domain); PythonEnvironment::gil acquire; - PYBIND11_OVERLOAD_PURE(double, ScalarField, getValue, pos.x(), pos.y(), pos.z()); } + + Vec3d getGradient(Vec3d& pos, int& domain) override { + PythonEnvironment::gil acquire; + + PYBIND11_OVERLOAD(Vec3d, ScalarField, getGradient, pos); + } + + void getHessian(Vec3d &pos, Mat3x3& h) override { + PythonEnvironment::gil acquire; + + PYBIND11_OVERLOAD(void, ScalarField, getHessian, pos, h); + } + }; void moduleAddScalarField(py::module &m) { @@ -70,11 +83,14 @@ void moduleAddScalarField(py::module &m) { py::cast(args[0]) + "'."); } } + //BindingBase::SetAttr(cc, key, value); } return ff; })); m.def("getValue", &ScalarField_Trampoline::getValue); + m.def("getGradient", &ScalarField_Trampoline::getGradient); + m.def("getHessian", &ScalarField_Trampoline::getHessian); } }