From 597bf61b468987b945617be0820cbe052464be4d Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Tue, 15 Sep 2026 14:36:19 +0200 Subject: [PATCH 1/8] Clean DiscreteGridField by moving the MHD loader's code into a separate file. --- CMakeLists.txt | 2 + src/SofaImplicitField/MHD.cpp | 162 +++++++++++++++ src/SofaImplicitField/MHD.h | 34 +++ .../components/geometry/DiscreteGridField.cpp | 195 ++---------------- .../components/geometry/DiscreteGridField.h | 50 ++--- .../initSofaImplicitField.cpp | 4 +- 6 files changed, 236 insertions(+), 211 deletions(-) create mode 100644 src/SofaImplicitField/MHD.cpp create mode 100644 src/SofaImplicitField/MHD.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 02f508f..8ab0a82 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,6 +13,7 @@ set(HEADER_FILES ${SOFAIMPLICITFIELD_SRC}/config.h.in ${SOFAIMPLICITFIELD_SRC}/initSofaImplicitField.h ${SOFAIMPLICITFIELD_SRC}/MarchingCube.h + ${SOFAIMPLICITFIELD_SRC}/MHD.h # This is backward compatibility ${SOFAIMPLICITFIELD_SRC}/deprecated/SphereSurface.h @@ -32,6 +33,7 @@ set(HEADER_FILES set(SOURCE_FILES ${SOFAIMPLICITFIELD_SRC}/initSofaImplicitField.cpp ${SOFAIMPLICITFIELD_SRC}/MarchingCube.cpp + ${SOFAIMPLICITFIELD_SRC}/MHD.cpp ## This is a backward compatibility.. ${SOFAIMPLICITFIELD_SRC}/deprecated/SphereSurface.cpp diff --git a/src/SofaImplicitField/MHD.cpp b/src/SofaImplicitField/MHD.cpp new file mode 100644 index 0000000..8d74e25 --- /dev/null +++ b/src/SofaImplicitField/MHD.cpp @@ -0,0 +1,162 @@ +/****************************************************************************** +* 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::loader +{ +using sofa::type::Vec3d; +using sofa::type::Vec3u; + +bool loadGridFromMHD( const std::string& filename_, Vec3d& imgMin, Vec3d& spacing, Vec3u& imgSize, float*& imgData) +{ + const char* filename = filename_.c_str(); + imgMin={0.0, 0.0, 0.0}; + imgSize={0, 0, 0}; + spacing={0, 0, 0}; + + char buffer[1024]; + char *value; + bool dataFileSpecified = false; + float f0, f1, f2; + int i0, i1, i2; + char dataFile[1024]; + + // read header file + std::ifstream header( filename ); + if (!header.is_open()) return false; + while (!header.eof()) + { + header.getline( buffer, 1024 ); + if (strncmp( buffer, "ObjectType", 10 ) == 0) + { + value = strchr( buffer, '=' )+1; while (*value==' ') value++; + if (strncmp( value, "Image", 5 ) != 0) + { + printf( "ERROR: Object is no image.\n" ); + return false; + } + } + else if (strncmp( buffer, "NDims", 5 ) == 0) + { + value = strchr( buffer, '=' )+1; while (*value==' ') value++; + if (*value != '3') + { + printf( "ERROR: Wrong number of dimensions.\n" ); + return false; + } + } + else if (strncmp( buffer, "BinaryData ", 11 ) == 0 || strncmp( buffer, "BinaryData=", 11 ) == 0) + { + value = strchr( buffer, '=' )+1; while (*value==' ') value++; + if (strncmp( value, "True", 4 ) != 0) + { + printf( "ERROR: Data is not binary.\n" ); + return false; + } + } + else if (strncmp( buffer, "CompressedData", 14 ) == 0) + { + value = strchr( buffer, '=' )+1; while (*value==' ') value++; + if (strncmp( value, "False", 5 ) != 0) + { + printf( "ERROR: Data is compressed.\n" ); + return false; + } + } + else if (strncmp( buffer, "TransformMatrix", 15 ) == 0) + { + value = strchr( buffer, '=' )+1; while (*value==' ') value++; + if (strncmp( value, "1 0 0 0 1 0 0 0 1", 17 ) != 0) + { + printf( "ERROR: Unsupported transform matrix.\n" ); + return false; + } + } + else if (strncmp( buffer, "Offset", 6 ) == 0) + { + value = strchr( buffer, '=' )+1; while (*value==' ') value++; + sscanf( value, "%f %f %f", &f0, &f1, &f2 ); + imgMin[0]=f0; imgMin[1]=f1; imgMin[2]=f2; + printf( "Image offset = %f %f %f\n", imgMin[0], imgMin[1], imgMin[2] ); + } + else if (strncmp( buffer, "ElementSpacing", 14 ) == 0) + { + value = strchr( buffer, '=' )+1; while (*value==' ') value++; + sscanf( value, "%f %f %f", &f0, &f1, &f2 ); + spacing[0]=f0; spacing[1]=f1; spacing[2]=f2; + printf( "Image spacing = %f %f %f\n", spacing[0], spacing[1], spacing[2] ); + } + else if (strncmp( buffer, "DimSize", 7 ) == 0) + { + value = strchr( buffer, '=' )+1; while (*value==' ') value++; + sscanf( value, "%d %d %d", &i0, &i1, &i2 ); + imgSize[0]=i0; imgSize[1]=i1; imgSize[2]=i2; + printf( "Image size = %i %i %i\n", imgSize[0], imgSize[1], imgSize[2] ); + } + else if (strncmp( buffer, "ElementType", 11 ) == 0) + { + value = strchr( buffer, '=' )+1; while (*value==' ') value++; + if (strncmp( value, "MET_FLOAT", 9) != 0) + { + printf( "ERROR: Datatype is not supported.\n" ); + return false; + } + } + } + header.close(); + + // read data file + if (!dataFileSpecified) + { + // change extension to .raw + strncpy( dataFile, filename, sizeof(dataFile) - 1 ); + dataFile[sizeof(dataFile) - 1] = '\0'; + size_t lenWithoutExt = strlen( filename ); + if (lenWithoutExt >= 3) + lenWithoutExt -= 3; + if (lenWithoutExt < sizeof(dataFile) - 4) + { + dataFile[lenWithoutExt] = '\0'; + strncat( dataFile, "raw", sizeof(dataFile) - lenWithoutExt - 1 ); + } + else + { + printf( "Warning: filename too long to replace extension, keeping '%s'\n", dataFile ); + } + } + std::ifstream data( dataFile, std::ios_base::binary|std::ios_base::in ); + if (!data.is_open()) return false; + unsigned int numVoxels = imgSize[0]* imgSize[1]* imgSize[2]; + imgData = new float[numVoxels]; + data.read( (char*)imgData, numVoxels*sizeof(float) ); + if (data.bad()) return false; + data.close(); + return false; +} + + +} diff --git a/src/SofaImplicitField/MHD.h b/src/SofaImplicitField/MHD.h new file mode 100644 index 0000000..f2777a3 --- /dev/null +++ b/src/SofaImplicitField/MHD.h @@ -0,0 +1,34 @@ +/****************************************************************************** +* 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 + +namespace sofaimplicitfield::loader +{ +using sofa::type::Vec3d; +using sofa::type::Vec3u; + +bool loadGridFromMHD(const std::string& filename_, Vec3d& imgMin, Vec3d& imgSize, Vec3u& spacing, float*& imgData); +} + diff --git a/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp b/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp index 83327db..db3d33e 100644 --- a/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp +++ b/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp @@ -19,63 +19,26 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ -#include #include -#include -using sofa::core::RegisterObject ; - #include +#include +#include +using sofa::core::RegisterObject ; -namespace sofa::component::geometry::_discretegrid_ +namespace sofa::component::geometry { -/** -DiscreteGridField::DiscreteGridField() - : in_filename(initData(&in_filename,"filename","filename")) - , in_nx(initData(&in_nx,0,"nx","in_nx")) - , in_ny(initData(&in_ny,0,"ny","in_ny")) - , in_nz(initData(&in_nz,0,"nz","in_nz")) - , in_scale(initData(&in_scale,0.0,"scale","in_scale")) - , in_sampling(initData(&in_sampling,0.0,"sampling","in_sampling")) -{ -} - - - -void DiscreteGridField::init() -{ - if(in_nx.getValue()==0 && in_nz.getValue()==0 && in_nz.getValue()==0) { - d_componentState.setValue(ComponentState::Invalid); - msg_error() << "uninitialized grid"; - } - else if(in_filename.isSet() == false) { - d_componentState.setValue(ComponentState::Invalid) - msg_error() << "unset filename"; - } - else { - pmin.set(0,0,-5.0); - pmax.set(27,27,5.0); - loadGrid(in_scale.getValue(),in_sampling.getValue(),in_nx.getValue(),in_ny.getValue(),in_nz.getValue(),pmin,pmax); - } - - d_componentState.setValue(ComponentState::Valid) -} -*/ - DiscreteGridField::DiscreteGridField() : ScalarField(), d_distanceMapHeader( initData( &d_distanceMapHeader, "file", "MHD file for the distance map" ) ), d_maxDomains( initData( &d_maxDomains, 1, "maxDomains", "Number of domains available for caching" ) ), - dx( initData( &dx, 0.0, "dx", "x translation" ) ), - dy( initData( &dy, 0.0, "dy", "y translation" ) ), - dz( initData( &dz, 0.0, "dz", "z translation" ) ) + d_position(initData( &d_position, {0.0,0.0,0.0}, "position", "The position in world space of the grid" ) ) { m_usedDomains = 0; m_imgData = nullptr; } - DiscreteGridField::~DiscreteGridField() { if (m_imgData) @@ -85,14 +48,12 @@ DiscreteGridField::~DiscreteGridField() } } - ///used to set a name in tests void DiscreteGridField::setFilename(const std::string& name) { d_distanceMapHeader.setValue(name); } - void DiscreteGridField::init() { m_domainCache.resize( d_maxDomains.getValue() ); @@ -100,110 +61,12 @@ void DiscreteGridField::init() if (ok) printf( "Successfully loaded distance map.\n" ); } - bool DiscreteGridField::loadGridFromMHD( const char *filename ) { - m_imgMin[0]=m_imgMin[1]=m_imgMin[2] = 0; - m_spacing[0]=m_spacing[1]=m_spacing[2] = 1; - m_imgSize[0]=m_imgSize[1]=m_imgSize[2] = 0; + bool loadSucceeded = sofaimplicitfield::loader::loadGridFromMHD(filename, m_imgMin, m_spacing, m_imgSize, m_imgData); - char buffer[1024]; - char *value; - bool dataFileSpecified = false; - float f0, f1, f2; - int i0, i1, i2; - char dataFile[1024]; - - // read header file - std::ifstream header( filename ); - if (!header.is_open()) return false; - while (!header.eof()) - { - header.getline( buffer, 1024 ); - if (strncmp( buffer, "ObjectType", 10 ) == 0) - { - value = strchr( buffer, '=' )+1; while (*value==' ') value++; - if (strncmp( value, "Image", 5 ) != 0) - { - printf( "ERROR: Object is no image.\n" ); - return false; - } - } - else if (strncmp( buffer, "NDims", 5 ) == 0) - { - value = strchr( buffer, '=' )+1; while (*value==' ') value++; - if (*value != '3') - { - printf( "ERROR: Wrong number of dimensions.\n" ); - return false; - } - } - else if (strncmp( buffer, "BinaryData ", 11 ) == 0 || strncmp( buffer, "BinaryData=", 11 ) == 0) - { - value = strchr( buffer, '=' )+1; while (*value==' ') value++; - if (strncmp( value, "True", 4 ) != 0) - { - printf( "ERROR: Data is not binary.\n" ); - return false; - } - } - else if (strncmp( buffer, "CompressedData", 14 ) == 0) - { - value = strchr( buffer, '=' )+1; while (*value==' ') value++; - if (strncmp( value, "False", 5 ) != 0) - { - printf( "ERROR: Data is compressed.\n" ); - return false; - } - } - else if (strncmp( buffer, "TransformMatrix", 15 ) == 0) - { - value = strchr( buffer, '=' )+1; while (*value==' ') value++; - if (strncmp( value, "1 0 0 0 1 0 0 0 1", 17 ) != 0) - { - printf( "ERROR: Unsupported transform matrix.\n" ); - return false; - } - } - else if (strncmp( buffer, "Offset", 6 ) == 0) - { - value = strchr( buffer, '=' )+1; while (*value==' ') value++; - sscanf( value, "%f %f %f", &f0, &f1, &f2 ); - m_imgMin[0]=f0; m_imgMin[1]=f1; m_imgMin[2]=f2; - printf( "Image offset = %f %f %f\n", m_imgMin[0], m_imgMin[1], m_imgMin[2] ); - } - else if (strncmp( buffer, "ElementSpacing", 14 ) == 0) - { - value = strchr( buffer, '=' )+1; while (*value==' ') value++; - sscanf( value, "%f %f %f", &f0, &f1, &f2 ); - m_spacing[0]=f0; m_spacing[1]=f1; m_spacing[2]=f2; - printf( "Image spacing = %f %f %f\n", m_spacing[0], m_spacing[1], m_spacing[2] ); - } - else if (strncmp( buffer, "DimSize", 7 ) == 0) - { - value = strchr( buffer, '=' )+1; while (*value==' ') value++; - sscanf( value, "%d %d %d", &i0, &i1, &i2 ); - m_imgSize[0]=i0; m_imgSize[1]=i1; m_imgSize[2]=i2; - printf( "Image size = %i %i %i\n", m_imgSize[0], m_imgSize[1], m_imgSize[2] ); - } - else if (strncmp( buffer, "ElementType", 11 ) == 0) - { - value = strchr( buffer, '=' )+1; while (*value==' ') value++; - if (strncmp( value, "MET_FLOAT", 9) != 0) - { - printf( "ERROR: Datatype is not supported.\n" ); - return false; - } - } - /* Don't allow variable names for problems with correct file paths! - else if (strncmp( buffer, "ElementDataFile", 11 ) == 0) { - value = strchr( buffer, '=' )+1; while (*value==' ') value++; - strncpy( dataFile, value, sizeof(dataFile) - 1 ); - dataFile[sizeof(dataFile) - 1] = '\0'; - dataFileSpecified = true; - }*/ - } - header.close(); + if(!loadSucceeded) + return false; // init remaining variables for (int d=0; d<3; d++) @@ -221,43 +84,20 @@ bool DiscreteGridField::loadGridFromMHD( const char *filename ) m_deltaOfs[6] = m_deltaOfs[2] + sliceSize; m_deltaOfs[7] = m_deltaOfs[3] + sliceSize; - // read data file - if (!dataFileSpecified) - { - // change extension to .raw - strncpy( dataFile, filename, sizeof(dataFile) - 1 ); - dataFile[sizeof(dataFile) - 1] = '\0'; - size_t lenWithoutExt = strlen( filename ); - if (lenWithoutExt >= 3) - lenWithoutExt -= 3; - if (lenWithoutExt < sizeof(dataFile) - 4) - { - dataFile[lenWithoutExt] = '\0'; - strncat( dataFile, "raw", sizeof(dataFile) - lenWithoutExt - 1 ); - } - else - { - printf( "Warning: filename too long to replace extension, keeping '%s'\n", dataFile ); - } - } - std::ifstream data( dataFile, std::ios_base::binary|std::ios_base::in ); - if (!data.is_open()) return false; - unsigned int numVoxels = m_imgSize[0]*m_imgSize[1]*m_imgSize[2]; - m_imgData = new float[numVoxels]; - data.read( (char*)m_imgData, numVoxels*sizeof(float) ); - if (data.bad()) return false; - data.close(); return true; } void DiscreteGridField::updateCache( DomainCache *cache, Vec3d& pos ) { cache->insideImg = true; - for (int d=0; d<3; d++) if (pos[d]=m_imgMax[d]) + for (int d=0; d<3; d++) + { + if (pos[d]=m_imgMax[d]) { cache->insideImg = false; break; } + } if (cache->insideImg) { int voxMinPos[3]; @@ -303,7 +143,6 @@ void DiscreteGridField::updateCache( DomainCache *cache, Vec3d& pos ) } } - int DiscreteGridField::getNextDomain() { // while we have free domains always return the next one, afterwards always use the last one @@ -311,14 +150,11 @@ int DiscreteGridField::getNextDomain() return m_usedDomains-1; } - double DiscreteGridField::getValue( Vec3d &transformedPos, int &domain ) { // use translation - Vec3d pos; - pos[0] = transformedPos[0] - dx.getValue(); - pos[1] = transformedPos[1] - dy.getValue(); - pos[2] = transformedPos[2] - dz.getValue(); + Vec3d pos = d_position.getValue(); + // find cache domain and check if it needs an update DomainCache *cache; if (domain < 0) @@ -359,7 +195,6 @@ double DiscreteGridField::getValue( Vec3d &transformedPos, int &domain ) return res; } - double DiscreteGridField::getValue( Vec3d &transformedPos ) { static int domain=-1; @@ -373,4 +208,4 @@ void registerDiscreteGridField(sofa::core::ObjectFactory* factory) .add< DiscreteGridField >()); } -} ///namespace sofa::component::geometry::_discretegrid_ +} \ No newline at end of file diff --git a/src/SofaImplicitField/components/geometry/DiscreteGridField.h b/src/SofaImplicitField/components/geometry/DiscreteGridField.h index a98e897..4ada320 100644 --- a/src/SofaImplicitField/components/geometry/DiscreteGridField.h +++ b/src/SofaImplicitField/components/geometry/DiscreteGridField.h @@ -19,28 +19,23 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ -#ifndef SOFAIMPLICITFIELD_COMPONENT_DISCRETEGRIDFIELD_H -#define SOFAIMPLICITFIELD_COMPONENT_DISCRETEGRIDFIELD_H +#pragma once #include #include #include -namespace sofa +namespace sofa::component::geometry { -namespace component +namespace { - -namespace geometry -{ - -namespace _discretegrid_ -{ - +using sofa::type::Vec3; +using sofa::type::Vec3u; using sofa::type::Vec3d; +} -class SOFA_SOFAIMPLICITFIELD_API DomainCache +class SOFA_SOFAIMPLICITFIELD_API DomainCache { public: bool insideImg; // shows if the domain lies inside the valid image region or outside @@ -72,28 +67,25 @@ class SOFA_SOFAIMPLICITFIELD_API DiscreteGridField : public virtual ScalarField sofa::core::objectmodel::DataFileName d_distanceMapHeader; Data< int > d_maxDomains; ///< Number of domains available for caching - Data< double > dx; ///< x translation - Data< double > dy; ///< y translation - Data< double > dz; ///< z translation - int m_usedDomains; // number of domains already given out - unsigned int m_imgSize[3]; // number of voxels - double m_spacing[3]; // physical distance between two neighboring voxels - double m_scale[3]; // (1/spacing) - double m_imgMin[3], m_imgMax[3]; // physical locations of the centers of both corner voxels + //Data< double > dx; ///< x translation + //Data< double > dy; ///< y translation + //Data< double > dz; ///< z translation + + Data< Vec3 > d_position; + + Vec3u m_imgSize; // number of voxels + Vec3d m_spacing; // physical distance between two neighboring voxels + Vec3d m_scale; // (1/spacing) + Vec3d m_imgMin; + Vec3d m_imgMax; // physical locations of the centers of both corner voxels float *m_imgData; // raw data + + int m_usedDomains; // number of domains already given out unsigned int m_deltaOfs[8]; // offsets to define 8 corners of cube for interpolation std::vector m_domainCache; }; -} /// namespace _discretegrid_ -using _discretegrid_::DiscreteGridField ; - -} /// namespace geometry - -} /// namespace component - -} /// namespace sofa +} -#endif diff --git a/src/SofaImplicitField/initSofaImplicitField.cpp b/src/SofaImplicitField/initSofaImplicitField.cpp index 9e3a2d5..f7da52a 100644 --- a/src/SofaImplicitField/initSofaImplicitField.cpp +++ b/src/SofaImplicitField/initSofaImplicitField.cpp @@ -46,7 +46,7 @@ namespace sofa::component::container { extern void registerInterpolatedImplicitSurface(sofa::core::ObjectFactory* factory); } -namespace sofa::component::geometry::_discretegrid_ +namespace sofa::component::geometry { extern void registerDiscreteGridField(sofa::core::ObjectFactory* factory); } @@ -108,7 +108,7 @@ void registerObjects(sofa::core::ObjectFactory* factory) sofa::component::geometry::_StarShapedField_::registerStarShapedField(factory); sofaimplicitfield::mapping::registerImplicitSurfaceMapping(factory); sofa::component::container::registerInterpolatedImplicitSurface(factory); - sofa::component::geometry::_discretegrid_::registerDiscreteGridField(factory); + sofa::component::geometry::registerDiscreteGridField(factory); sofaimplicitfield::component::engine::registerFieldToSurfaceMesh(factory); } From 3bf4a0e95bf46a49898dfa039c0ce66d05b2b452 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Tue, 22 Sep 2026 16:27:51 +0200 Subject: [PATCH 2/8] First working version of a modernized DiscreteGridField. The sampling is now done in component called GridSampler while the DiscreteGridField is there just to do the tri linear interpolation. --- CMakeLists.txt | 2 + .../example-grid-generation-from-implicit.py | 50 +++ examples/python/xshape/primitives.py | 24 +- python/src/Binding_ScalarField.cpp | 5 + src/SofaImplicitField/MarchingCube.cpp | 5 +- .../components/engine/GridSampler.cpp | 180 ++++++++ .../components/engine/GridSampler.h | 72 +++ .../components/geometry/DiscreteGridField.cpp | 410 +++++++++++++----- .../components/geometry/DiscreteGridField.h | 79 ++-- .../initSofaImplicitField.cpp | 6 +- 10 files changed, 674 insertions(+), 159 deletions(-) create mode 100644 examples/python/example-grid-generation-from-implicit.py create mode 100644 src/SofaImplicitField/components/engine/GridSampler.cpp create mode 100644 src/SofaImplicitField/components/engine/GridSampler.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 8ab0a82..6a7d436 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -21,6 +21,7 @@ set(HEADER_FILES ${SOFAIMPLICITFIELD_SRC}/deprecated/InterpolatedImplicitSurface.h # This is a backward compatibility file toward DiscreteGridField ${SOFAIMPLICITFIELD_SRC}/components/engine/FieldToSurfaceMesh.h + ${SOFAIMPLICITFIELD_SRC}/components/engine/GridSampler.h ${SOFAIMPLICITFIELD_SRC}/components/geometry/BottleField.h ${SOFAIMPLICITFIELD_SRC}/components/geometry/DiscreteGridField.h ${SOFAIMPLICITFIELD_SRC}/components/geometry/SphericalField.h @@ -40,6 +41,7 @@ set(SOURCE_FILES ${SOFAIMPLICITFIELD_SRC}/deprecated/InterpolatedImplicitSurface.cpp ${SOFAIMPLICITFIELD_SRC}/components/engine/FieldToSurfaceMesh.cpp + ${SOFAIMPLICITFIELD_SRC}/components/engine/GridSampler.cpp ${SOFAIMPLICITFIELD_SRC}/components/geometry/BottleField.cpp ${SOFAIMPLICITFIELD_SRC}/components/geometry/ScalarField.cpp ${SOFAIMPLICITFIELD_SRC}/components/geometry/DiscreteGridField.cpp diff --git a/examples/python/example-grid-generation-from-implicit.py b/examples/python/example-grid-generation-from-implicit.py new file mode 100644 index 0000000..55d9cc0 --- /dev/null +++ b/examples/python/example-grid-generation-from-implicit.py @@ -0,0 +1,50 @@ +import Sofa +import SofaImplicitField + +from Sofa.Types import RGBAColor +from SofaTypes.SofaTypes import Vec3d +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 grids from two different scalar field and visualize them as mesh + using mesh extraction from implicit field. + """ + root.addObject("RequiredPlugin", pluginName="SofaImplicitField") + + root.addObject(DrawController()) + + ########################### The same field is used for both grids ################## + if True: + f1 = root.addObject( + Difference(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)) + ) + else: + #f1 = root.addObject(SpatialField(name="sphere", axis=0)) + f1 = root.addObject(Sphere(name="sphere", center=[0,0,0],radius=1.0)) + + ########################### GridGeneration ################## + g = root.addChild("Grid_10x10x10") + m1 = g.addObject("GridSampler", name="sampler", min=[-2,-2,-2], max=[2,2,2], resolution=[255,255,255]) + m1.field.setLinkedBase(f1) + + e1 = g.addObject("DiscreteGridField", name="grid1", buffer=m1.buffer.linkpath) + + m2 = g.addObject("FieldToSurfaceMesh", name="polygonizer1", + field=e1.linkpath, min=[-2,-2,-2], max=[2,2,2], step=0.01) + r = g.addObject("OglModel", name="renderer", + position=g.polygonizer1.points.linkpath, + triangles=g.polygonizer1.triangles.linkpath) + diff --git a/examples/python/xshape/primitives.py b/examples/python/xshape/primitives.py index b4d6422..e828740 100644 --- a/examples/python/xshape/primitives.py +++ b/examples/python/xshape/primitives.py @@ -18,14 +18,6 @@ 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 - def getValues(self, positions, results): """This version of the overrides the getValues so that we fetch the data once""" results[:] = numpy.linalg.norm(positions - self.center.value, axis=1) - self.radius.value @@ -54,4 +46,20 @@ def getValues(self, positions, results): outside = numpy.linalg.norm(numpy.maximum(q, 0.0), axis=1) inside = numpy.minimum(numpy.max(q, axis=1), 0.0) results[:] = outside + inside - r + return results + +class SpatialField(ScalarField): + def __init__(self, *args, **kwargs): + ScalarField.__init__(self, *args, **kwargs) + + self.addData("axis", type="int",value=kwargs.get("axis", 0), default=0, help="axis of the spatial field among {0, 1, 2}", group="Geometry") + + def getValue(self, pos): + """This version is of very low performance as there are a huge amount of call to the python side""" + x,y,z = pos + return x + + def getValues(self, positions, results): + """This version of the overrides the getValues so that we fetch the data once""" + results[:] = positions[:, self.axis.value] return results \ No newline at end of file diff --git a/python/src/Binding_ScalarField.cpp b/python/src/Binding_ScalarField.cpp index e9ae6d3..f3266d3 100644 --- a/python/src/Binding_ScalarField.cpp +++ b/python/src/Binding_ScalarField.cpp @@ -210,6 +210,11 @@ void moduleAddScalarField(py::module &m) { self->getHessian(pos, result); return result; }); + + /// register the PointSetTopologyModifier binding in the downcasting subsystem + PythonFactory::registerType([](sofa::core::objectmodel::Base* object) { + return py::cast(dynamic_cast(object)); + }); } } diff --git a/src/SofaImplicitField/MarchingCube.cpp b/src/SofaImplicitField/MarchingCube.cpp index 4a7d009..42245de 100644 --- a/src/SofaImplicitField/MarchingCube.cpp +++ b/src/SofaImplicitField/MarchingCube.cpp @@ -34,7 +34,7 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, const Vec3d& gridmin, const Vec3d& gridmax, std::function&, std::vector&)> getFieldValueAt, SeqCoord& tmpPoints, SeqTriangles& tmpTriangles) -{ +{ int nx = floor((gridmax.x() - gridmin.x()) * invStep) + 1 ; int ny = floor((gridmax.y() - gridmin.y()) * invStep) + 1 ; int nz = floor((gridmax.z() - gridmin.z()) * invStep) + 1 ; @@ -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; @@ -61,8 +60,6 @@ void MarchingCube::generateSurfaceMesh(const double isoval, const double mstep, const int dx = 1; const int dy = nx; - z = 0; - 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) diff --git a/src/SofaImplicitField/components/engine/GridSampler.cpp b/src/SofaImplicitField/components/engine/GridSampler.cpp new file mode 100644 index 0000000..1f4f16d --- /dev/null +++ b/src/SofaImplicitField/components/engine/GridSampler.cpp @@ -0,0 +1,180 @@ +/****************************************************************************** +* 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 +using sofa::core::RegisterObject; +using sofa::core::visual::VisualParams; + +namespace sofaimplicitfield::component::engine +{ + +GridSampler::GridSampler() + : d_resolution(initData(&d_resolution, Vec3u(10, 10, 10), "resolution", "Number of samples in each dimension")) + , d_min(initData(&d_min, Vec3d(-1.0, -1.0, -1.0), "min", "Minimum corner of the sampling grid")) + , d_max(initData(&d_max, Vec3d(1.0, 1.0, 1.0), "max", "Maximum corner of the sampling grid")) + , d_buffer(initData(&d_buffer, "buffer", "The data buffer olding the values")) + , l_field(initLink("field", "The scalar field to sample")) +{ + addUpdateCallback("sample", {&d_resolution, &d_min, &d_max}, [this](const sofa::core::DataTracker&) + { + Vec3u resolution = d_resolution.getValue(); + Vec3d min = d_min.getValue(); + Vec3d max = d_max.getValue(); + + updateInternalBuffer(resolution, min, max); + sampleField(); + return core::objectmodel::ComponentState::Valid; + }, {}); +} + +GridSampler::~GridSampler() +{ +} + +void GridSampler::init() +{ + if (!l_field.get()) + { + msg_error() << "Missing scalar field to sample"; + d_componentState = core::objectmodel::ComponentState::Invalid; + return; + } + + updateInternalBuffer(d_resolution.getValue(), d_min.getValue(), d_max.getValue()); + sampleField(); + d_componentState = core::objectmodel::ComponentState::Valid; +} + +void GridSampler::notifyLinkSet(BaseLink* link, Base* base) +{ + if(link==&l_field){ + std::cout << "THE FIELD HAS HCNAGD" << std::endl; + } +} + +void GridSampler::updateInternalBuffer(const Vec3u& resolution, const Vec3d& min, const Vec3d& max) +{ + auto buffer = sofa::helper::getWriteOnlyAccessor(d_buffer); + buffer->resize(min, max, resolution); +} + +void GridSampler::sampleField() +{ + auto field = l_field.get(); + if (!field) + { + msg_error() << "No scalar field linked"; + return; + } + + auto buffer = sofa::helper::getWriteAccessor(d_buffer); + auto resolution = buffer->resolution; + auto min = buffer->min; + auto max = buffer->max; + + if (resolution[0] <= 0 || resolution[1] <= 0 || resolution[2] <= 0) + { + msg_error() << "Resolution must be greater than 0 in all dimensions"; + return; + } + + if (min[0] >= max[0] || min[1] >= max[1] || min[2] >= max[2]) + { + msg_error() << "min must be less than max in all dimensions"; + return; + } + + std::cout << getPathName() << " sampling the grid " << std::endl; + auto data = buffer->data; + auto spacing = buffer->spacing; + + std::cout << " : " << buffer->data << std::endl; + + int rx = resolution.x(); + int rxy = resolution.y() * rx; + + std::vector positions; + std::vector results; + positions.reserve(rxy); + results.reserve(rxy); + + for (unsigned int z = 0; z < resolution[2]; z++) + { + positions.clear(); + results.clear(); + for (unsigned int y = 0; y < resolution[1]; y++) + { + for (unsigned int x = 0; x < resolution[0]; x++) + { + positions.emplace_back(min + Vec3d( + static_cast(x) * spacing[0], + static_cast(y) * spacing[1], + static_cast(z) * spacing[2] + )); + } + } + field->getValues(positions, results); + for (unsigned int y = 0; y < resolution[1]; y++) + { + for (unsigned int x = 0; x < resolution[0]; x++) + { + auto index=[rxy,rx](int x,int y, int z){ return z * rxy + y * rx + x; }; + auto rindex=[rx](int x,int y){ return y * rx + x; }; + + data[index(x,y,z)] = results[rindex(x,y)]; + } + } + } + std::cout << "SAMPLING DONE "<< std::endl; +} + +void GridSampler::computeBBox(const core::ExecParams* params, bool onlyVisible) +{ + if (onlyVisible) return; + + Vec3d min = d_min.getValue(); + Vec3d max = d_max.getValue(); + + sofa::core::objectmodel::BaseComponent::computeBBox(params, onlyVisible); + + f_bbox.setValue({min, max}); +} + +void GridSampler::draw(const VisualParams* vparams) +{ + if (!vparams || !vparams->displayFlags().getShowBehaviorModels()) + return; + + if (isComponentStateInvalid()) + return; +} + +void registerGridSampler(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Samples a ScalarField on a 3D grid and stores the result in a DiscreteGridField.") + .add< GridSampler >()); +} + +} diff --git a/src/SofaImplicitField/components/engine/GridSampler.h b/src/SofaImplicitField/components/engine/GridSampler.h new file mode 100644 index 0000000..1fcf5a1 --- /dev/null +++ b/src/SofaImplicitField/components/engine/GridSampler.h @@ -0,0 +1,72 @@ +/****************************************************************************** +* 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 + +//////////////////////////////////////////////////////////////////////////////////////////////////// +namespace sofaimplicitfield::component::engine +{ + +namespace{ + using namespace sofa; + using sofa::core::objectmodel::BaseComponent; + using sofa::core::visual::VisualParams; + using sofa::type::Vec3d; + using sofa::type::Vec3u; + using sofa::component::geometry::ScalarField; + using sofa::component::geometry::DiscreteGridField; +} + +class GridSampler : public BaseComponent +{ +public: + SOFA_CLASS(GridSampler, BaseComponent); + + void init() override; + void draw(const VisualParams* params) override; + + Data d_resolution; + Data d_min; + Data d_max; + + Data d_buffer; + + void notifyLinkSet(BaseLink*, Base*) override; + +protected: + SingleLink l_field; + +protected: + GridSampler(); + virtual ~GridSampler(); + +private: + void computeBBox(const core::ExecParams* params, bool onlyVisible = false) override; + void updateInternalBuffer(const Vec3u& resolution, const Vec3d& min, const Vec3d& max); + void sampleField(); +}; + +} + diff --git a/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp b/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp index db3d33e..9696e54 100644 --- a/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp +++ b/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp @@ -21,184 +21,343 @@ ******************************************************************************/ #include #include +#include #include +#include +#include + #include using sofa::core::RegisterObject ; +#include + namespace sofa::component::geometry { DiscreteGridField::DiscreteGridField() : ScalarField(), d_distanceMapHeader( initData( &d_distanceMapHeader, "file", "MHD file for the distance map" ) ), - d_maxDomains( initData( &d_maxDomains, 1, "maxDomains", "Number of domains available for caching" ) ), - d_position(initData( &d_position, {0.0,0.0,0.0}, "position", "The position in world space of the grid" ) ) + d_min(initData( &d_min, {-0.5,-0.5,-0.5}, "min", "The min positions in world space" ) ), + d_max(initData( &d_max, {0.5 ,0.5 ,0.5}, "max", "The max positions in world space" ) ), + d_resolution(initData( &d_resolution, {10,10,10}, "resolution", "The resolution of each axis of the grid" ) ), + d_buffer(initData(&d_buffer, "buffer", "The data buffer olding the values")), + d_debugDraw(initData(&d_debugDraw, false, "debugDraw", "show the values on the grid")) { - m_usedDomains = 0; - m_imgData = nullptr; + addUpdateCallback("updateFromRMM",{&d_resolution, &d_min, &d_max},[this](const sofa::core::DataTracker&){ + + /// Update the internal buffers when d_resolution change + std::cout << "WE ARE GOING TO UPDATE FROM DATA " << d_resolution.getCounter() << std::endl; + + auto buffer = sofa::helper::getWriteOnlyAccessor(d_buffer); + buffer->resize(d_min.getValue(), d_max.getValue(), d_resolution.getValue()); + internalUpdate(buffer); + refreshTrackers(); + + return sofa::core::objectmodel::ComponentState::Valid; + }, {}); + + addUpdateCallback("updateFromData",{&d_buffer},[this](const sofa::core::DataTracker&){ + + /// Update the internal buffers when d_resolution change + std::cout << "WE ARE GOING TO UPDATE FROM DATA BUFFER: " << d_buffer.getCounter() << std::endl; + + /// Update the internal state + auto buffer = sofa::helper::getReadAccessor(d_buffer); + internalUpdate(buffer); + + /// Propagate the change to the other inputs. + d_resolution.setValue(buffer->resolution); + d_min.setValue(buffer->min); + d_max.setValue(buffer->max); + + refreshTrackers(); + + return sofa::core::objectmodel::ComponentState::Valid; + }, {}); } DiscreteGridField::~DiscreteGridField() { - if (m_imgData) - { - delete[] m_imgData; - m_imgData = nullptr; - } } -///used to set a name in tests -void DiscreteGridField::setFilename(const std::string& name) +// Clean the tracker so we don't call the update mechanisme twice +void DiscreteGridField::refreshTrackers() { - d_distanceMapHeader.setValue(name); + for(auto& tracker : m_internalEngine){ + tracker.second.cleanDirty(); + } } void DiscreteGridField::init() { - m_domainCache.resize( d_maxDomains.getValue() ); - bool ok = loadGridFromMHD( d_distanceMapHeader.getFullPath().c_str() ); - if (ok) printf( "Successfully loaded distance map.\n" ); + if(d_distanceMapHeader.isSet()) + { + bool ok = loadGridFromMHD( d_distanceMapHeader.getFullPath().c_str() ); + if (ok){ + msg_info() << "Successfully loaded a distance map from file."; + d_componentState = sofa::core::objectmodel::ComponentState::Valid; + } + else{ + d_componentState = sofa::core::objectmodel::ComponentState::Invalid; + } + return; + } + std::cout << getPathName() << " init" << std::endl; + + if(!d_buffer.isSet()){ + internalResize(d_resolution.getValue(), d_min.getValue(), d_max.getValue()); + }else{ + const auto& buffer = sofa::helper::getReadAccessor(d_buffer); + std::cout << getPathName() << " buffer " << buffer->min << ", " << buffer->max << " and " << buffer->data << std::endl; + internalUpdate(buffer); + } + + std::cout << getPathName() << " init done" << std::endl; + d_componentState = sofa::core::objectmodel::ComponentState::Valid; } bool DiscreteGridField::loadGridFromMHD( const char *filename ) { - bool loadSucceeded = sofaimplicitfield::loader::loadGridFromMHD(filename, m_imgMin, m_spacing, m_imgSize, m_imgData); + auto buffer = sofa::helper::getWriteOnlyAccessor(d_buffer); + + bool loadSucceeded = sofaimplicitfield::loader::loadGridFromMHD(filename, + buffer->min, buffer->spacing, + buffer->resolution, buffer->data); if(!loadSucceeded) return false; // init remaining variables - for (int d=0; d<3; d++) - { - m_scale[d] = 1.0/m_spacing[d]; - m_imgMax[d] = m_imgMin[d] + (double)(m_imgSize[d]-1)*m_spacing[d]; - } - m_deltaOfs[0] = 0; - m_deltaOfs[1] = 1; - m_deltaOfs[2] = m_imgSize[0]; - m_deltaOfs[3] = m_imgSize[0]+1; - unsigned int sliceSize = m_imgSize[0]*m_imgSize[1]; - m_deltaOfs[4] = m_deltaOfs[0] + sliceSize; - m_deltaOfs[5] = m_deltaOfs[1] + sliceSize; - m_deltaOfs[6] = m_deltaOfs[2] + sliceSize; - m_deltaOfs[7] = m_deltaOfs[3] + sliceSize; + // for (int d=0; d<3; d++) + // { + // scaling[d] = 1.0/spacing[d]; + // max[d] = min[d] + (double)(resolution[d]-1)*spacing[d]; + // } + // m_deltaOfs[0] = 0; + // m_deltaOfs[1] = 1; + // m_deltaOfs[2] = resolution[0]; + // m_deltaOfs[3] = resolution[0]+1; + + // unsigned int sliceSize = resolution[0]*resolution[1]; + // m_deltaOfs[4] = m_deltaOfs[0] + sliceSize; + // m_deltaOfs[5] = m_deltaOfs[1] + sliceSize; + // m_deltaOfs[6] = m_deltaOfs[2] + sliceSize; + // m_deltaOfs[7] = m_deltaOfs[3] + sliceSize; return true; } -void DiscreteGridField::updateCache( DomainCache *cache, Vec3d& pos ) +void MemoryBuffer::resize(const Vec3d& min_, const Vec3d& max_, const Vec3u& resolution_) { - cache->insideImg = true; - for (int d=0; d<3; d++) - { - if (pos[d]=m_imgMax[d]) - { - cache->insideImg = false; - break; - } - } - if (cache->insideImg) - { - int voxMinPos[3]; - for (int d=0; d<3; d++) - { - voxMinPos[d] = (int)(m_scale[d] * (pos[d]-m_imgMin[d])); - cache->bbMin[d] = m_spacing[d]*(double)voxMinPos[d] + m_imgMin[d]; - cache->bbMax[d] = cache->bbMin[d] + m_spacing[d]; - } - unsigned int ofs = voxMinPos[0] + m_imgSize[0]*(voxMinPos[1] + m_imgSize[1]*voxMinPos[2]); - cache->val[0] = m_imgData[ofs]; - for (int i=1; i<8; i++) cache->val[i] = m_imgData[ofs+m_deltaOfs[i]]; - } - else - { - // init bounding box to be as large as possible to prevent unnecessary cache updates while outside image - const double MIN=-10e6, MAX=10e6; - int voxMappedPos[3]; - for (int d=0; d<3; d++) + resolution = resolution_; + min = min_; + max = max_; + + auto newImgSize = resolution.x() * resolution.y() * resolution.z(); + if(size != newImgSize){ + std::cout << "Resizeing memory buffer " << data << " to " << newImgSize << std::endl; + if (data) { - if (pos[d] < m_imgMin[d]) - { - cache->bbMin[d] = MIN; - cache->bbMax[d] = m_imgMin[d]; - voxMappedPos[d] = 0; - } - else if (pos[d] >= m_imgMax[d]) - { - cache->bbMin[d] = m_imgMax[d]; - cache->bbMax[d] = MAX; - voxMappedPos[d] = m_imgSize[d]-1; - } - else - { - cache->bbMin[d] = MIN; - cache->bbMax[d] = MAX; - voxMappedPos[d] = (int)(m_scale[d] * (pos[d]-m_imgMin[d])); - } + delete[] data; } - unsigned int ofs = voxMappedPos[0] + m_imgSize[0]*(voxMappedPos[1] + m_imgSize[1]*voxMappedPos[2]); - // if cache lies outside image, the returned distance is not updated anymore, instead this boundary value is returned - cache->val[0] = m_imgData[ofs] + m_spacing[0]+m_spacing[1]+m_spacing[2]; + data = new float[newImgSize]; + size = newImgSize; } + std::cout << "Resizeing done... " << data << std::endl; + + spacing = (max - min).linearDivision(Vec3d{ + static_cast(resolution[0] - 1 > 0 ? resolution[0] - 1 : 1), + static_cast(resolution[1] - 1 > 0 ? resolution[1] - 1 : 1), + static_cast(resolution[2] - 1 > 0 ? resolution[2] - 1 : 1) + }); + + scaling = Vec3d{1.0,1.0,1.0}.linearDivision(spacing); } -int DiscreteGridField::getNextDomain() +void DiscreteGridField::internalUpdate(const MemoryBuffer& buffer) { - // while we have free domains always return the next one, afterwards always use the last one - if (m_usedDomains < (int)m_domainCache.size()) m_usedDomains++; - return m_usedDomains-1; + const auto& resolution = buffer.resolution; + + unsigned int sliceSize = resolution[0] * resolution[1]; + m_deltaOfs[0] = 0; + m_deltaOfs[1] = 1; + m_deltaOfs[2] = resolution[0]; + m_deltaOfs[3] = resolution[0] + 1; + m_deltaOfs[4] = sliceSize; + m_deltaOfs[5] = sliceSize + 1; + m_deltaOfs[6] = sliceSize + resolution[0]; + m_deltaOfs[7] = sliceSize + resolution[0] + 1; + + d_resolution.setValue(resolution); + d_min.setValue(buffer.min); + d_max.setValue(buffer.max); } -double DiscreteGridField::getValue( Vec3d &transformedPos, int &domain ) +void DiscreteGridField::internalResize(const Vec3u& resolution, const Vec3d& min, const Vec3d& max) { - // use translation - Vec3d pos = d_position.getValue(); + auto buffer = sofa::helper::getWriteOnlyAccessor(d_buffer); + buffer->resize(min,max,resolution); + internalUpdate(buffer); +} - // find cache domain and check if it needs an update - DomainCache *cache; - if (domain < 0) - { - domain = getNextDomain(); - cache = &(m_domainCache[domain]); - updateCache( cache, pos ); - } - else +bool DiscreteGridField::empty() +{ + auto buffer = sofa::helper::getReadAccessor(d_buffer); + return buffer->data == nullptr; +} + +void DiscreteGridField::draw(const sofa::core::visual::VisualParams* params) +{ + if(!d_debugDraw.getValue()) + return; + + auto dt = params->drawTool(); + const auto& buffer = d_buffer.getValue(); + auto& resolution = buffer.resolution; + auto& min = buffer.min; + auto& data = buffer.data; + auto& spacing = buffer.spacing; + + auto index = [resolution](unsigned int x, unsigned int y, unsigned int z) { return x + resolution.x() * y + (resolution.x() * resolution.y() * z); }; + + for(unsigned int x=0;xbbMin[d] || pos[d]>cache->bbMax[d]) + for(unsigned int z=0;zdraw3DText(pos, 0.1, type::RGBAColor::cyan(), s.str().c_str()); } } } +} - // if cache lies outside image, the returned distance is not updated anymore, instead this boundary value is returned - if (!cache->insideImg) return cache->val[0]; +void DiscreteGridField::getValues(const std::vector& positions, std::vector& results) +{ + auto buffer = sofa::helper::getReadAccessor(d_buffer); + if(buffer->data==nullptr) + return; - // use trilinear interpolation on cached cube - double weight[3]; - for (int d=0; d<3; d++) - { - weight[d] = m_scale[d] * (pos[d]-cache->bbMin[d]); + const auto& resolution = buffer->resolution; + const auto& min = buffer->min; + const auto& scaling = buffer->scaling; + const auto& spacing = buffer->spacing; + const auto& data = buffer->data; + + auto getValue = [&min, &scaling, &resolution, &data](const Vec3d position){ + Vec3d localPosition = (position-min); + + Vec3d t = localPosition.linearProduct(scaling); + + // Compute the indices of voxels surrounding the position + unsigned int x0 = std::min((unsigned int)std::floor(t.x()), resolution.x() - 2 ); + unsigned int y0 = std::min((unsigned int)std::floor(t.y()), resolution.y() - 2 ); + unsigned int z0 = std::min((unsigned int)std::floor(t.z()), resolution.z() - 2 ); + + unsigned int x1 = x0 + 1; + unsigned int y1 = y0 + 1; + unsigned int z1 = z0 + 1; + + t = t-Vec3d{x0,y0,z0}; + + // Helper function to access the index + auto index = [resolution](unsigned int x, unsigned int y, unsigned int z) { return x + resolution.x() * y + (resolution.x() * resolution.y() * z); }; + + // Get the height raw values surrounding the position + const double c000 = data[index(x0, y0, z0)]; + const double c100 = data[index(x1, y0, z0)]; + const double c010 = data[index(x0, y1, z0)]; + const double c110 = data[index(x1, y1, z0)]; + const double c001 = data[index(x0, y0, z1)]; + const double c101 = data[index(x1, y0, z1)]; + const double c011 = data[index(x0, y1, z1)]; + const double c111 = data[index(x1, y1, z1)]; + + // X + const double c00 = c000 * (1.0 - t.x()) + c100 * t.x(); + const double c10 = c010 * (1.0 - t.x()) + c110 * t.x(); + const double c01 = c001 * (1.0 - t.x()) + c101 * t.x(); + const double c11 = c011 * (1.0 - t.x()) + c111 * t.x(); + + // Y + const double c0 = c00 * (1.0 - t.y()) + c10 * t.y(); + const double c1 = c01 * (1.0 - t.y()) + c11 * t.y(); + + // Z + const double value = c0 * (1.0 - t.z()) + c1 * t.z(); + + return value; + }; + + for(unsigned int i=0;ival[0]*a + cache->val[1]*b + cache->val[2]*c + cache->val[3]*d ) * (1.0-weight[2]) - + ( cache->val[4]*a + cache->val[5]*b + cache->val[6]*c + cache->val[7]*d ) * weight[2]; - - return res; + } -double DiscreteGridField::getValue( Vec3d &transformedPos ) +double DiscreteGridField::getValue(Vec3d &position, int &) { - static int domain=-1; - return getValue( transformedPos, domain ); + auto buffer = sofa::helper::getReadAccessor(d_buffer); + if(buffer->data==nullptr) + return -1.0; + + const auto& resolution = buffer->resolution; + const auto& min = buffer->min; + const auto& scaling = buffer->scaling; + const auto& spacing = buffer->spacing; + const auto& data = buffer->data; + + Vec3d localPosition = (position-min); + + // use trilinear interpolation to get the value at any location + Vec3d t = localPosition.linearProduct(scaling); + + // Compute the indices of voxels surrounding the position + unsigned int x0 = std::min((unsigned int)std::floor(t.x()), resolution.x() - 2 ); + unsigned int y0 = std::min((unsigned int)std::floor(t.y()), resolution.y() - 2 ); + unsigned int z0 = std::min((unsigned int)std::floor(t.z()), resolution.z() - 2 ); + + unsigned int x1 = x0 + 1; + unsigned int y1 = y0 + 1; + unsigned int z1 = z0 + 1; + + t = t-Vec3d{x0,y0,z0}; + + // Helper function to access the index + auto index = [resolution](unsigned int x, unsigned int y, unsigned int z) { return x + resolution.x() * y + (resolution.x() * resolution.y() * z); }; + + // Get the height raw values surrounding the position + const double c000 = data[index(x0, y0, z0)]; + const double c100 = data[index(x1, y0, z0)]; + const double c010 = data[index(x0, y1, z0)]; + const double c110 = data[index(x1, y1, z0)]; + const double c001 = data[index(x0, y0, z1)]; + const double c101 = data[index(x1, y0, z1)]; + const double c011 = data[index(x0, y1, z1)]; + const double c111 = data[index(x1, y1, z1)]; + + // X + const double c00 = c000 * (1.0 - t.x()) + c100 * t.x(); + const double c10 = c010 * (1.0 - t.x()) + c110 * t.x(); + const double c01 = c001 * (1.0 - t.x()) + c101 * t.x(); + const double c11 = c011 * (1.0 - t.x()) + c111 * t.x(); + + // Y + const double c0 = c00 * (1.0 - t.y()) + c10 * t.y(); + const double c1 = c01 * (1.0 - t.y()) + c11 * t.y(); + + // Z + const double value = c0 * (1.0 - t.z()) + c1 * t.z(); + + //std::cout << "POSITION " << position << "("<< x0 << "," << y0 << "," << z0 << " and "<< t << ") value " << c000 << std::endl; + return value; } // Register in the Factory @@ -208,4 +367,17 @@ void registerDiscreteGridField(sofa::core::ObjectFactory* factory) .add< DiscreteGridField >()); } +} + +namespace sofa::core::objectmodel +{ +/// Specialization for MemoryBuffer +template<> bool Data::read( const std::string&) { return false; } +template<> void Data::printValue( std::ostream& ) const {} +template<> std::string Data::getValueString() const { + std::stringstream tmp; + auto& buffer = getValue(); + tmp << "Grid<" << buffer.resolution.x() << "," << buffer.resolution.y() << "," << buffer.resolution.z() << "> @"< std::string Data::getDefaultValueString() const { return ""; } } \ No newline at end of file diff --git a/src/SofaImplicitField/components/geometry/DiscreteGridField.h b/src/SofaImplicitField/components/geometry/DiscreteGridField.h index 4ada320..0a4ccf9 100644 --- a/src/SofaImplicitField/components/geometry/DiscreteGridField.h +++ b/src/SofaImplicitField/components/geometry/DiscreteGridField.h @@ -25,6 +25,7 @@ #include #include + namespace sofa::component::geometry { @@ -35,57 +36,83 @@ using sofa::type::Vec3u; using sofa::type::Vec3d; } -class SOFA_SOFAIMPLICITFIELD_API DomainCache +class MemoryBuffer { public: - bool insideImg; // shows if the domain lies inside the valid image region or outside - Vec3d bbMin, bbMax; // bounding box (min and max) of the domain - double val[8]; // corner values of the domain + Vec3d spacing; + Vec3d scaling; + Vec3d min; + Vec3d max; + Vec3u resolution; + unsigned int size{0}; + float* data{nullptr}; + + void resize(const Vec3d& min_, const Vec3d& max_, const Vec3u& resolution_); }; + class SOFA_SOFAIMPLICITFIELD_API DiscreteGridField : public virtual ScalarField { - public: SOFA_CLASS(DiscreteGridField, ScalarField); -public: DiscreteGridField(); ~DiscreteGridField() override; void init() override; + void draw(const sofa::core::visual::VisualParams*) override; - virtual double getValue( Vec3d &transformedPos ); - double getValue( Vec3d &transformedPos, int &domain ) override; - int getDomain( Vec3d &pos, int ref_domain ) override { (void)pos; return ref_domain; } + double getValue(Vec3d& position, int& domain) override; + void getValues(const std::vector& positions, std::vector& results) override; - void setFilename(const std::string& filename) ; bool loadGridFromMHD( const char *filename ) ; - void updateCache( DomainCache *cache, Vec3d& pos ); - int getNextDomain(); - sofa::core::objectmodel::DataFileName d_distanceMapHeader; - Data< int > d_maxDomains; ///< Number of domains available for caching - //Data< double > dx; ///< x translation - //Data< double > dy; ///< y translation - //Data< double > dz; ///< z translation + Data d_min; // bounding box (min) + Data d_max; // bounding box (max) + Data d_resolution; // resolution of the grid along each axis - Data< Vec3 > d_position; + Data d_buffer; - Vec3u m_imgSize; // number of voxels - Vec3d m_spacing; // physical distance between two neighboring voxels - Vec3d m_scale; // (1/spacing) - Vec3d m_imgMin; - Vec3d m_imgMax; // physical locations of the centers of both corner voxels - float *m_imgData; // raw data + Data d_debugDraw; - int m_usedDomains; // number of domains already given out unsigned int m_deltaOfs[8]; // offsets to define 8 corners of cube for interpolation - std::vector m_domainCache; + bool empty(); + + class Modifier + { + public: + Modifier(DiscreteGridField* field){self=field;} + ~Modifier(){ self->getComponentState(); } + Modifier& resize(const Vec3u& resolution, + const Vec3d& gridMin, const Vec3d& gridMax){ + + self->d_min.setValue(gridMin); + self->d_max.setValue(gridMax); + self->d_resolution.setValue(resolution); + return *this; } //< resize the grid and re-allocate the buffers + private: + DiscreteGridField* self; + }; + Modifier modify(){ return Modifier(this); } + void refreshTrackers(); + +private: + void internalUpdate(const MemoryBuffer& buffer); + void internalResize(const Vec3u& resolution, const Vec3d& min, const Vec3d& max); + }; } +namespace sofa::core::objectmodel +{ + +/// Specialization for MemoryBuffer +template<> bool Data::read( const std::string&); +template<> void Data::printValue( std::ostream& ) const; +template<> std::string Data::getValueString() const; +template<> std::string Data::getDefaultValueString() const; +} diff --git a/src/SofaImplicitField/initSofaImplicitField.cpp b/src/SofaImplicitField/initSofaImplicitField.cpp index f7da52a..974d8ff 100644 --- a/src/SofaImplicitField/initSofaImplicitField.cpp +++ b/src/SofaImplicitField/initSofaImplicitField.cpp @@ -14,7 +14,7 @@ * * * 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 * @@ -52,7 +52,8 @@ namespace sofa::component::geometry } namespace sofaimplicitfield::component::engine { -extern void registerFieldToSurfaceMesh(sofa::core::ObjectFactory* factory); + extern void registerFieldToSurfaceMesh(sofa::core::ObjectFactory* factory); + extern void registerGridSampler(sofa::core::ObjectFactory* factory); } namespace sofaimplicitfield @@ -110,6 +111,7 @@ void registerObjects(sofa::core::ObjectFactory* factory) sofa::component::container::registerInterpolatedImplicitSurface(factory); sofa::component::geometry::registerDiscreteGridField(factory); sofaimplicitfield::component::engine::registerFieldToSurfaceMesh(factory); + sofaimplicitfield::component::engine::registerGridSampler(factory); } } /// sofaimplicitfield From 7a8a8b3887fd312fcae81deabe7bbdd312eb4214 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 23 Sep 2026 13:25:45 +0200 Subject: [PATCH 3/8] Change signature of ScalarField::getValue to have const for position. --- python/src/Binding_ScalarField.cpp | 6 +- .../components/geometry/BottleField.cpp | 58 ++++++++++--------- .../components/geometry/BottleField.h | 6 +- .../components/geometry/DiscreteGridField.cpp | 6 +- .../components/geometry/DiscreteGridField.h | 2 +- .../components/geometry/ScalarField.cpp | 9 +-- .../components/geometry/ScalarField.h | 22 +++---- .../components/geometry/SphericalField.cpp | 13 +++-- .../components/geometry/SphericalField.h | 6 +- .../components/geometry/StarShapedField.cpp | 6 +- .../components/geometry/StarShapedField.h | 6 +- 11 files changed, 71 insertions(+), 69 deletions(-) diff --git a/python/src/Binding_ScalarField.cpp b/python/src/Binding_ScalarField.cpp index f3266d3..b61030b 100644 --- a/python/src/Binding_ScalarField.cpp +++ b/python/src/Binding_ScalarField.cpp @@ -103,7 +103,7 @@ class ScalarField_Trampoline : public ScalarField { return py::str(py::cast(this).get_type().attr("__name__")); } - double getValue(Vec3& pos, int& domain) override + double getValue(const Vec3& pos, int& domain) override { SOFA_UNUSED(domain); PythonEnvironment::gil acquire; @@ -129,7 +129,7 @@ class ScalarField_Trampoline : public ScalarField { auto o = override(vector_to_numpy(positions), vector_to_numpy(results)); } - Vec3 getGradient(Vec3& pos, int& domain) override + Vec3 getGradient(const Vec3& pos, int& domain) override { SOFA_UNUSED(domain); PythonEnvironment::gil acquire; @@ -137,7 +137,7 @@ class ScalarField_Trampoline : public ScalarField { PYBIND11_OVERLOAD(Vec3, ScalarField, getGradient, pos); } - void getHessian(Vec3 &pos, Mat3x3& h) override + void getHessian(const Vec3 &pos, Mat3x3& h) override { /// The implementation is a bit more complex compared to getGradient. This is because we change de signature between the c++ API and the python one. PythonEnvironment::gil acquire; diff --git a/src/SofaImplicitField/components/geometry/BottleField.cpp b/src/SofaImplicitField/components/geometry/BottleField.cpp index 11b9bb1..32fa51f 100644 --- a/src/SofaImplicitField/components/geometry/BottleField.cpp +++ b/src/SofaImplicitField/components/geometry/BottleField.cpp @@ -76,11 +76,12 @@ double BottleField::innerLength(Vec3d& Pos) m_excentricity*(Pos[2] - m_center[2])*(Pos[2] - m_center[2])); } -double BottleField::getValue(Vec3d& Pos, int& domain) +double BottleField::getValue(const Vec3d& pos_, int& domain) { SOFA_UNUSED(domain) ; - double resultSphereOuter = this->outerLength(Pos) - m_radius ; - double resultEllipsoidInner = this->innerLength(Pos) - m_ellipsoidRadius; + Vec3d pos = pos_; + double resultSphereOuter = this->outerLength(pos) - m_radius ; + double resultEllipsoidInner = this->innerLength(pos) - m_ellipsoidRadius; double result = std::max(resultSphereOuter,-resultEllipsoidInner); @@ -90,24 +91,24 @@ double BottleField::getValue(Vec3d& Pos, int& domain) return result; } -Vec3d BottleField::getGradient(Vec3d &Pos, int &domain) +Vec3d BottleField::getGradient(const Vec3d& pos_, int &domain) { SOFA_UNUSED(domain); Vec3d g; - - double LsphereOuter = this->outerLength(Pos) ; - double LEllipsoidInner = this->innerLength(Pos) ; + Vec3d pos = pos_; + double LsphereOuter = this->outerLength(pos) ; + double LEllipsoidInner = this->innerLength(pos) ; if (LsphereOuter - m_radius > - (LEllipsoidInner- m_ellipsoidRadius)){ - g[0] = (Pos[0] - m_center[0])/LsphereOuter; - g[1] = (Pos[1] - m_center[1])/LsphereOuter; - g[2] = (Pos[2] - m_center[2])/LsphereOuter; + g[0] = (pos[0] - m_center[0])/LsphereOuter; + g[1] = (pos[1] - m_center[1])/LsphereOuter; + g[2] = (pos[2] - m_center[2])/LsphereOuter; } else { - g[0] = -m_excentricity*(Pos[0] - m_center[0])/LEllipsoidInner; - g[1] = -(Pos[1] - (m_center[1]+m_shift))/LEllipsoidInner; - g[2] = -m_excentricity*(Pos[2] - m_center[2])/LEllipsoidInner; + g[0] = -m_excentricity*(pos[0] - m_center[0])/LEllipsoidInner; + g[1] = -(pos[1] - (m_center[1]+m_shift))/LEllipsoidInner; + g[2] = -m_excentricity*(pos[2] - m_center[2])/LEllipsoidInner; } @@ -121,34 +122,35 @@ Vec3d BottleField::getGradient(Vec3d &Pos, int &domain) return g; } -void BottleField::getHessian(Vec3d &Pos, Mat3x3& h) +void BottleField::getHessian(const Vec3d &pos_, Mat3x3& h) { - double LsphereOuter = this->outerLength(Pos) ; - double LEllipsoidInner = this->innerLength(Pos) ; + Vec3d pos = pos_; + double LsphereOuter = this->outerLength(pos) ; + double LEllipsoidInner = this->innerLength(pos) ; if (LsphereOuter - m_radius > - (LEllipsoidInner- m_ellipsoidRadius)) { double LsphereOuterSquare = LsphereOuter*LsphereOuter; double LsphereOuterCube = LsphereOuter*LsphereOuter*LsphereOuter; - h[0][0] = ( LsphereOuter - (Pos[0] - m_center[0])*(Pos[0] - m_center[0])/LsphereOuter )/LsphereOuterSquare ; - h[1][1] = ( LsphereOuter - (Pos[1] - m_center[1])*(Pos[1] - m_center[1])/LsphereOuter )/LsphereOuterSquare ; - h[2][2] = ( LsphereOuter - (Pos[2] - m_center[2])*(Pos[2] - m_center[2])/LsphereOuter )/LsphereOuterSquare ; + h[0][0] = ( LsphereOuter - (pos[0] - m_center[0])*(pos[0] - m_center[0])/LsphereOuter )/LsphereOuterSquare ; + h[1][1] = ( LsphereOuter - (pos[1] - m_center[1])*(pos[1] - m_center[1])/LsphereOuter )/LsphereOuterSquare ; + h[2][2] = ( LsphereOuter - (pos[2] - m_center[2])*(pos[2] - m_center[2])/LsphereOuter )/LsphereOuterSquare ; - h[0][1] = h[1][0] = - (Pos[0] - m_center[0])*(Pos[1] - m_center[1]) / LsphereOuterCube; - h[0][2] = h[2][0] = - (Pos[0] - m_center[0])*(Pos[2] - m_center[2]) / LsphereOuterCube; - h[1][2] = h[2][1] = - (Pos[2] - m_center[2])*(Pos[1] - m_center[1]) / LsphereOuterCube; + h[0][1] = h[1][0] = - (pos[0] - m_center[0])*(pos[1] - m_center[1]) / LsphereOuterCube; + h[0][2] = h[2][0] = - (pos[0] - m_center[0])*(pos[2] - m_center[2]) / LsphereOuterCube; + h[1][2] = h[2][1] = - (pos[2] - m_center[2])*(pos[1] - m_center[1]) / LsphereOuterCube; } else { double LEllipsoidInnerSquare = LEllipsoidInner*LEllipsoidInner; double LEllipsoidInnerCube = LEllipsoidInner*LEllipsoidInner*LEllipsoidInner; - h[0][0] = -m_excentricity*(LEllipsoidInner - m_excentricity*(Pos[0] - m_center[0])*(Pos[0] - m_center[0])/LEllipsoidInner )/LEllipsoidInnerSquare ; - h[1][1] = -(LEllipsoidInner - (Pos[1] - (m_center[1]+m_shift))*(Pos[1] - (m_center[1]+m_shift))/LEllipsoidInner )/LEllipsoidInnerSquare ; - h[2][2] = -m_excentricity*(LEllipsoidInner - m_excentricity*(Pos[2] - m_center[2])*(Pos[2] - m_center[2])/LEllipsoidInner )/LEllipsoidInnerSquare ; + h[0][0] = -m_excentricity*(LEllipsoidInner - m_excentricity*(pos[0] - m_center[0])*(pos[0] - m_center[0])/LEllipsoidInner )/LEllipsoidInnerSquare ; + h[1][1] = -(LEllipsoidInner - (pos[1] - (m_center[1]+m_shift))*(pos[1] - (m_center[1]+m_shift))/LEllipsoidInner )/LEllipsoidInnerSquare ; + h[2][2] = -m_excentricity*(LEllipsoidInner - m_excentricity*(pos[2] - m_center[2])*(pos[2] - m_center[2])/LEllipsoidInner )/LEllipsoidInnerSquare ; - h[0][1] = h[1][0] = m_excentricity*(Pos[0] - m_center[0])*(Pos[1] - (m_center[1]+m_shift)) / LEllipsoidInnerCube; - h[0][2] = h[2][0] = m_excentricity*m_excentricity*(Pos[0] - m_center[0])*(Pos[2] - m_center[2]) / LEllipsoidInnerCube; - h[1][2] = h[2][1] = m_excentricity*(Pos[2] - m_center[2])*(Pos[1] - (m_center[1]+m_shift)) / LEllipsoidInnerCube; + h[0][1] = h[1][0] = m_excentricity*(pos[0] - m_center[0])*(pos[1] - (m_center[1]+m_shift)) / LEllipsoidInnerCube; + h[0][2] = h[2][0] = m_excentricity*m_excentricity*(pos[0] - m_center[0])*(pos[2] - m_center[2]) / LEllipsoidInnerCube; + h[1][2] = h[2][1] = m_excentricity*(pos[2] - m_center[2])*(pos[1] - (m_center[1]+m_shift)) / LEllipsoidInnerCube; } if (m_inside) diff --git a/src/SofaImplicitField/components/geometry/BottleField.h b/src/SofaImplicitField/components/geometry/BottleField.h index 01249e4..f6c9fb3 100644 --- a/src/SofaImplicitField/components/geometry/BottleField.h +++ b/src/SofaImplicitField/components/geometry/BottleField.h @@ -50,9 +50,9 @@ class SOFA_SOFAIMPLICITFIELD_API BottleField : public ScalarField void reinit() override ; /// Inherited from ScalarField. - double getValue(Vec3d& Pos, int &domain) override ; - Vec3d getGradient(Vec3d &Pos, int& domain) override ; - void getHessian(Vec3d &Pos, Mat3x3& h) override; + double getValue(const Vec3d& Pos, int &domain) override ; + Vec3d getGradient(const Vec3d &Pos, int& domain) override ; + void getHessian(const Vec3d &Pos, Mat3x3& h) override; double outerLength(Vec3d& Pos); double innerLength(Vec3d& Pos); diff --git a/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp b/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp index 9696e54..9046430 100644 --- a/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp +++ b/src/SofaImplicitField/components/geometry/DiscreteGridField.cpp @@ -59,7 +59,7 @@ DiscreteGridField::DiscreteGridField() addUpdateCallback("updateFromData",{&d_buffer},[this](const sofa::core::DataTracker&){ - /// Update the internal buffers when d_resolution change + /// Update the internal buffers when d_buffer change std::cout << "WE ARE GOING TO UPDATE FROM DATA BUFFER: " << d_buffer.getCounter() << std::endl; /// Update the internal state @@ -248,7 +248,6 @@ void DiscreteGridField::getValues(const std::vector& positions, std::vect const auto& resolution = buffer->resolution; const auto& min = buffer->min; const auto& scaling = buffer->scaling; - const auto& spacing = buffer->spacing; const auto& data = buffer->data; auto getValue = [&min, &scaling, &resolution, &data](const Vec3d position){ @@ -302,7 +301,7 @@ void DiscreteGridField::getValues(const std::vector& positions, std::vect } -double DiscreteGridField::getValue(Vec3d &position, int &) +double DiscreteGridField::getValue(const Vec3d& position, int &) { auto buffer = sofa::helper::getReadAccessor(d_buffer); if(buffer->data==nullptr) @@ -356,7 +355,6 @@ double DiscreteGridField::getValue(Vec3d &position, int &) // Z const double value = c0 * (1.0 - t.z()) + c1 * t.z(); - //std::cout << "POSITION " << position << "("<< x0 << "," << y0 << "," << z0 << " and "<< t << ") value " << c000 << std::endl; return value; } diff --git a/src/SofaImplicitField/components/geometry/DiscreteGridField.h b/src/SofaImplicitField/components/geometry/DiscreteGridField.h index 0a4ccf9..b7fb948 100644 --- a/src/SofaImplicitField/components/geometry/DiscreteGridField.h +++ b/src/SofaImplicitField/components/geometry/DiscreteGridField.h @@ -62,7 +62,7 @@ class SOFA_SOFAIMPLICITFIELD_API DiscreteGridField : public virtual ScalarField void init() override; void draw(const sofa::core::visual::VisualParams*) override; - double getValue(Vec3d& position, int& domain) override; + double getValue(const Vec3d& position, int& domain) override; void getValues(const std::vector& positions, std::vector& results) override; bool loadGridFromMHD( const char *filename ) ; diff --git a/src/SofaImplicitField/components/geometry/ScalarField.cpp b/src/SofaImplicitField/components/geometry/ScalarField.cpp index de0396a..e90cb6b 100644 --- a/src/SofaImplicitField/components/geometry/ScalarField.cpp +++ b/src/SofaImplicitField/components/geometry/ScalarField.cpp @@ -46,9 +46,10 @@ void ScalarField::init() d_componentState.setValue(core::objectmodel::ComponentState::Valid); } -Vec3d ScalarField::getGradientByFinitDifference(Vec3d& pos, int& i) +Vec3d ScalarField::getGradientByFinitDifference(const Vec3d& pos_, int& i) { Vec3d Result; + Vec3d pos = pos_; double epsilon = d_epsilon.getValue(); pos[0] += epsilon; Result[0] = getValue(pos, i); @@ -78,12 +79,12 @@ void ScalarField::getValues(const std::vector& positions, 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); - inline Vec3d getGradient(Vec3d& pos) {int domain=-1; return getGradient(pos,domain); } - virtual void getHessian(Vec3d &Pos, Mat3x3& h); + virtual Vec3d getGradient(const Vec3d& pos, int& domain); + inline Vec3d getGradient(const Vec3d& pos) {int domain=-1; return getGradient(pos,domain); } + virtual void getHessian(const Vec3d &Pos, Mat3x3& h); /// Returns the value and the gradiant by evaluating one after an other. /// For some computation it is possible to implement more efficiently the computation /// By factoring the computing of the two...if you can do this please override this function. - virtual void getValueAndGradient(Vec3d& pos, double &value, Vec3d& grad, int& domain) ; - inline void getValueAndGradient(Vec3d& pos, double &value, Vec3d& grad) + virtual void getValueAndGradient(const Vec3d& pos, double &value, Vec3d& grad, int& domain) ; + inline void getValueAndGradient(const Vec3d& pos, double &value, Vec3d& grad) { int domain=-1; return getValueAndGradient(pos,value,grad,domain); diff --git a/src/SofaImplicitField/components/geometry/SphericalField.cpp b/src/SofaImplicitField/components/geometry/SphericalField.cpp index 9f492ab..87ed00e 100644 --- a/src/SofaImplicitField/components/geometry/SphericalField.cpp +++ b/src/SofaImplicitField/components/geometry/SphericalField.cpp @@ -47,12 +47,13 @@ void SphericalField::reinit() init(); } -double SphericalField::getValue(Vec3d& Pos, int& domain) +double SphericalField::getValue(const Vec3d& pos_, int& domain) { SOFA_UNUSED(domain) ; - double result = (Pos[0] - m_center[0])*(Pos[0] - m_center[0]) + - (Pos[1] - m_center[1])*(Pos[1] - m_center[1]) + - (Pos[2] - m_center[2])*(Pos[2] - m_center[2]) - + Vec3d pos = pos_; + double result = (pos[0] - m_center[0])*(pos[0] - m_center[0]) + + (pos[1] - m_center[1])*(pos[1] - m_center[1]) + + (pos[2] - m_center[2])*(pos[2] - m_center[2]) - m_radius * m_radius ; if(m_inside) result = -result; @@ -60,7 +61,7 @@ double SphericalField::getValue(Vec3d& Pos, int& domain) return result; } -Vec3d SphericalField::getGradient(Vec3d &Pos, int &domain) +Vec3d SphericalField::getGradient(const Vec3d &Pos, int &domain) { SOFA_UNUSED(domain); Vec3d g; @@ -80,7 +81,7 @@ Vec3d SphericalField::getGradient(Vec3d &Pos, int &domain) return g; } -void SphericalField::getValueAndGradient(Vec3d& Pos, double &value, Vec3d& /*grad*/, int& domain) +void SphericalField::getValueAndGradient(const Vec3d& Pos, double &value, Vec3d& /*grad*/, int& domain) { SOFA_UNUSED(domain); Vec3d g; diff --git a/src/SofaImplicitField/components/geometry/SphericalField.h b/src/SofaImplicitField/components/geometry/SphericalField.h index c77df18..87b1b91 100644 --- a/src/SofaImplicitField/components/geometry/SphericalField.h +++ b/src/SofaImplicitField/components/geometry/SphericalField.h @@ -52,9 +52,9 @@ class SOFA_SOFAIMPLICITFIELD_API SphericalField : public ScalarField void reinit() override ; /// Inherited from ScalarField. - double getValue(Vec3d& Pos, int &domain) override ; - Vec3d getGradient(Vec3d &Pos, int& domain) override ; - void getValueAndGradient(Vec3d& pos, double& val, Vec3d& grad, int& domain) override ; + double getValue(const Vec3d& Pos, int &domain) override ; + Vec3d getGradient(const Vec3d &Pos, int& domain) override ; + void getValueAndGradient(const Vec3d& pos, double& val, Vec3d& grad, int& domain) override ; using ScalarField::getValue ; using ScalarField::getGradient ; diff --git a/src/SofaImplicitField/components/geometry/StarShapedField.cpp b/src/SofaImplicitField/components/geometry/StarShapedField.cpp index 412c571..df09729 100644 --- a/src/SofaImplicitField/components/geometry/StarShapedField.cpp +++ b/src/SofaImplicitField/components/geometry/StarShapedField.cpp @@ -60,7 +60,7 @@ void StarShapedField::reinit() init(); } -double StarShapedField::getValue(Vec3d& Pos, int& domain) +double StarShapedField::getValue(const Vec3d& Pos, int& domain) { SOFA_UNUSED(domain) ; double length = sqrt((Pos[0] - m_center[0])*(Pos[0] - m_center[0]) + @@ -79,7 +79,7 @@ double StarShapedField::getValue(Vec3d& Pos, int& domain) return result; } -Vec3d StarShapedField::getGradient(Vec3d &Pos, int &domain) +Vec3d StarShapedField::getGradient(const Vec3d &Pos, int &domain) { SOFA_UNUSED(domain); Vec3d g; @@ -105,7 +105,7 @@ Vec3d StarShapedField::getGradient(Vec3d &Pos, int &domain) return g; } -void StarShapedField::getHessian(Vec3d &Pos, Mat3x3& h) +void StarShapedField::getHessian(const Vec3d &Pos, Mat3x3& h) { double length = sqrt((Pos[0] - m_center[0])*(Pos[0] - m_center[0]) + (Pos[1] - m_center[1])*(Pos[1] - m_center[1]) + diff --git a/src/SofaImplicitField/components/geometry/StarShapedField.h b/src/SofaImplicitField/components/geometry/StarShapedField.h index 3295bb6..90c5798 100644 --- a/src/SofaImplicitField/components/geometry/StarShapedField.h +++ b/src/SofaImplicitField/components/geometry/StarShapedField.h @@ -50,9 +50,9 @@ class SOFA_SOFAIMPLICITFIELD_API StarShapedField : public ScalarField void reinit() override ; /// Inherited from ScalarField. - double getValue(Vec3d& Pos, int &domain) override ; - Vec3d getGradient(Vec3d &Pos, int& domain) override ; - void getHessian(Vec3d &Pos, Mat3x3& h) override; + double getValue(const Vec3d& Pos, int &domain) override ; + Vec3d getGradient(const Vec3d &Pos, int& domain) override ; + void getHessian(const Vec3d &Pos, Mat3x3& h) override; using ScalarField::getValue ; using ScalarField::getGradient ; From 48d26f5414090830175dcc62459fd895f9dd3ea3 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Wed, 23 Sep 2026 13:26:07 +0200 Subject: [PATCH 4/8] Add a ScalarFieldMapping (WIP) --- CMakeLists.txt | 2 + examples/python/scalar-field-mapping.py | 43 +++++ .../components/mapping/ScalarFieldMapping.cpp | 160 ++++++++++++++++++ .../components/mapping/ScalarFieldMapping.h | 69 ++++++++ .../initSofaImplicitField.cpp | 2 + 5 files changed, 276 insertions(+) create mode 100644 examples/python/scalar-field-mapping.py create mode 100644 src/SofaImplicitField/components/mapping/ScalarFieldMapping.cpp create mode 100644 src/SofaImplicitField/components/mapping/ScalarFieldMapping.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 6a7d436..b1dce2a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,6 +29,7 @@ set(HEADER_FILES ${SOFAIMPLICITFIELD_SRC}/components/geometry/StarShapedField.h ${SOFAIMPLICITFIELD_SRC}/components/mapping/ImplicitSurfaceMapping.h ${SOFAIMPLICITFIELD_SRC}/components/mapping/ImplicitSurfaceMapping.inl + ${SOFAIMPLICITFIELD_SRC}/components/mapping/ScalarFieldMapping.h ) set(SOURCE_FILES @@ -48,6 +49,7 @@ set(SOURCE_FILES ${SOFAIMPLICITFIELD_SRC}/components/geometry/SphericalField.cpp ${SOFAIMPLICITFIELD_SRC}/components/geometry/StarShapedField.cpp ${SOFAIMPLICITFIELD_SRC}/components/mapping/ImplicitSurfaceMapping.cpp + ${SOFAIMPLICITFIELD_SRC}/components/mapping/ScalarFieldMapping.cpp ) set(EXTRA_FILES diff --git a/examples/python/scalar-field-mapping.py b/examples/python/scalar-field-mapping.py new file mode 100644 index 0000000..b8f2730 --- /dev/null +++ b/examples/python/scalar-field-mapping.py @@ -0,0 +1,43 @@ +import Sofa +import Sofa.Core +from SofaImplicitField import ScalarField +import numpy + +def Model(): + particules = Sofa.Core.Node("Particles") + particules.addObject("MechanicalObject", name="state", template="Vec3", position=[[-1,0,0],[0,0,0],[1,0,0]]) + particules.addObject("UniformMass", name="mass", totalMass=1.0) + + child = particules.addChild("Child") + child.addObject("MechanicalObject", name="state", template="Vec1", position=[0,0,0]) + child.addObject("SphericalField", name="field", center=[0,0,0], radius=1.0) + child.addObject("ScalarFieldMapping", name="mapping", + input=particules.state.linkpath, + output=child.state.linkpath, field=child.field.linkpath) + + return particules + +def createScene(root): + """In this scene we create a ScalarField of spherical shape and a mechanical object of 3D particles. + The field and the particules are used as input of a ScalarFieldMapping, mapping the particules to a 1D child space. + This child space is then used in a constraints. + """ + root.addObject('RequiredPlugin', pluginName='Sofa.Component.AnimationLoop') # Needed to use components [FreeMotionAnimationLoop] + root.addObject('RequiredPlugin', pluginName='Sofa.Component.Constraint.Lagrangian.Model') # Needed to use components [StopperLagrangianConstraint] + root.addObject('RequiredPlugin', pluginName='Sofa.Component.Constraint.Lagrangian.Solver') # Needed to use components [BlockGaussSeidelConstraintSolver] + root.addObject('RequiredPlugin', pluginName='Sofa.Component.StateContainer') # Needed to use components [MechanicalObject] + root.addObject("RequiredPlugin", pluginName="SofaImplicitField") + root.addObject('RequiredPlugin', pluginName='Sofa.Component.Mass') # Needed to use components [UniformMass] + + root.addObject("FreeMotionAnimationLoop") + root.addObject("BlockGaussSeidelConstraintSolver", name="solver", tolerance=1e-6, maxIterations=1000) + + root.addObject("EulerImplicitSolver", name="odesolver", rayleighStiffness=0.1, rayleighMass=0.1) + root.addObject("SparseLDLSolver", name="linearSolver", template="CompressedRowSparseMatrixd") + + model = root.addChild(Model()) + model.state.showObject = True + model.state.showObjectScale = 5.0 + + model.Child.addObject("StopperLagrangianConstraint", name="constraint", min=0.0, max=1.0, index=1) + model.Child.addObject("GenericConstraintCorrection", name="correction", linearSolver=root.linearSolver.linkpath, ODESolver=root.odesolver.linkpath) diff --git a/src/SofaImplicitField/components/mapping/ScalarFieldMapping.cpp b/src/SofaImplicitField/components/mapping/ScalarFieldMapping.cpp new file mode 100644 index 0000000..3ffd076 --- /dev/null +++ b/src/SofaImplicitField/components/mapping/ScalarFieldMapping.cpp @@ -0,0 +1,160 @@ +/****************************************************************************** +* 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 +using sofa::core::RegisterObject; + +namespace sofaimplicitfield::mapping +{ + +ScalarFieldMapping::ScalarFieldMapping() : + l_field(initLink("field", "The scalar field to sample")) +{ +} + +ScalarFieldMapping::~ScalarFieldMapping() +{ +} + +void ScalarFieldMapping::init() +{ + if(l_field.get() == nullptr){ + msg_error() << "The field is missing. Cannot work properly without it."; + d_componentState = core::objectmodel::ComponentState::Invalid; + return; + } + d_componentState = core::objectmodel::ComponentState::Valid; +} + +/// +void ScalarFieldMapping::apply( const MechanicalParams* mparams, OutDataVecCoord& out_, const InDataVecCoord& in_) +{ + SOFA_UNUSED(mparams); + if(!isComponentStateValid()) + return; + + auto in = sofa::helper::getReadAccessor(in_); + auto out = sofa::helper::getWriteOnlyAccessor(out_); + auto field = l_field.get(); + int domain{0}; + + for(Size i=0; igetValue(in[i],domain)); + } +} + +/// This method must be reimplemented by all mappings. +void ScalarFieldMapping::applyJ( const MechanicalParams* mparams, OutDataVecDeriv& dy_, const InDataVecDeriv& dx_) +{ + SOFA_UNUSED(mparams); + if(!isComponentStateValid()) + return; + auto x = fromModel->readPositions(); + auto dx = sofa::helper::getReadAccessor(dx_); + auto dy = sofa::helper::getWriteOnlyAccessor(dy_); + auto field = l_field.get(); + + for (Size i = 0; i < x.size(); ++i) + { + const Vec3d grad = field->getGradient(x[i]); + dy[i] = sofa::type::dot(grad, dx[i]); + } +} + +/// This method must be reimplemented by all mappings. +void ScalarFieldMapping::applyJT( const MechanicalParams* mparams, InDataVecDeriv& dx_, const OutDataVecDeriv& dy_) +{ + SOFA_UNUSED(mparams); + if(!isComponentStateValid()) + return; + auto x = fromModel->readPositions(); + auto dx = sofa::helper::getWriteOnlyAccessor(dx_); + auto dy = sofa::helper::getReadAccessor(dy_); + auto field = l_field.get(); + + for (Size i = 0; i < dx.size(); ++i) + { + const Vec3d grad = field->getGradient(x[i]); + dx[i] += grad * dy[i].x(); //< Because dy is in R so there is only one value + } +} + +void ScalarFieldMapping::applyJT( const ConstraintParams* mparams, InDataMatrixDeriv& dx_, const OutDataMatrixDeriv& dy_) +{ + SOFA_UNUSED(mparams); + if(!isComponentStateValid()) + return; + auto x = fromModel->readPositions(); + auto dx = sofa::helper::getWriteOnlyAccessor(dx_); + auto dy = sofa::helper::getReadAccessor(dy_); + auto field = l_field.get(); + + for (Size i = 0; i < dx.size(); ++i) + { + const Vec3d grad = field->getGradient(x[i]); + dx[i] += grad * dy[i].x(); //< Because dy is in R so there is only one value + } +} + + +void ScalarFieldMapping::buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) +{ + if(!isComponentStateValid()) + return; + + const auto childForces = this->toModel->readTotalForces(); + const auto dJdx = matrices->getMappingDerivativeIn(this->fromModel).withRespectToPositionsIn(this->fromModel); + const auto x = fromModel->readPositions(); + + auto field = l_field.get(); + + for (Size i = 0; i < x.size(); ++i) + { + const auto f = childForces[i]; + + type::Mat3x3d H; + field->getHessian(x[i], H); + + const sofa::type::Mat3x3d Ki = H * f.x(); // again les forces sont des "forces" 1d. + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + dJdx(3*i + a, 3*i + b) += Ki[a][b]; + } + } + } +} + +void registerScalarFieldMapping(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Maps a positional field to its scalar field values.") + .add< ScalarFieldMapping >()); +} + +} diff --git a/src/SofaImplicitField/components/mapping/ScalarFieldMapping.h b/src/SofaImplicitField/components/mapping/ScalarFieldMapping.h new file mode 100644 index 0000000..45cf983 --- /dev/null +++ b/src/SofaImplicitField/components/mapping/ScalarFieldMapping.h @@ -0,0 +1,69 @@ +/****************************************************************************** +* 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 + +//////////////////////////////////////////////////////////////////////////////////////////////////// +namespace sofaimplicitfield::mapping +{ + +namespace{ + using namespace sofa; + using sofa::core::objectmodel::BaseComponent; + using sofa::core::visual::VisualParams; + using sofa::core::ConstraintParams; + using sofa::type::Vec3d; + using sofa::type::Vec3u; + using sofa::component::geometry::ScalarField; + using sofa::core::Mapping; + using sofa::core::BaseMapping; + using sofa::core::MechanicalParams; + using sofa::core::MultiVecDerivId; + using sofa::core::ConstMultiVecDerivId; +} + +class ScalarFieldMapping : public Mapping +{ +public: + SOFA_CLASS(ScalarFieldMapping, + SOFA_TEMPLATE2(Mapping, sofa::defaulttype::Vec3Types, sofa::defaulttype::Vec1Types)); + + void init() override; + + ScalarFieldMapping(); + ~ScalarFieldMapping(); + + void apply( const MechanicalParams* mparams, OutDataVecCoord& out, const InDataVecCoord& in) override; + void applyJ( const MechanicalParams* mparams, OutDataVecDeriv& out, const InDataVecDeriv& in) override; + void applyJT( const MechanicalParams* mparams, InDataVecDeriv& out, const OutDataVecDeriv& in) override; + void applyJT( const ConstraintParams* /* mparams */, InDataMatrixDeriv& /* out */, const OutDataMatrixDeriv& /* in */) override; + + void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; + + SingleLink l_field; +}; + +} + diff --git a/src/SofaImplicitField/initSofaImplicitField.cpp b/src/SofaImplicitField/initSofaImplicitField.cpp index 974d8ff..590f1de 100644 --- a/src/SofaImplicitField/initSofaImplicitField.cpp +++ b/src/SofaImplicitField/initSofaImplicitField.cpp @@ -41,6 +41,7 @@ namespace sofa::component::geometry::_StarShapedField_ namespace sofaimplicitfield::mapping { extern void registerImplicitSurfaceMapping(sofa::core::ObjectFactory* factory); + extern void registerScalarFieldMapping(sofa::core::ObjectFactory* factory); } namespace sofa::component::container { @@ -108,6 +109,7 @@ void registerObjects(sofa::core::ObjectFactory* factory) sofa::component::geometry::_sphericalfield_::registerSphericalField(factory); sofa::component::geometry::_StarShapedField_::registerStarShapedField(factory); sofaimplicitfield::mapping::registerImplicitSurfaceMapping(factory); + sofaimplicitfield::mapping::registerScalarFieldMapping(factory); sofa::component::container::registerInterpolatedImplicitSurface(factory); sofa::component::geometry::registerDiscreteGridField(factory); sofaimplicitfield::component::engine::registerFieldToSurfaceMesh(factory); From 8200e90984c9d5b73e001c751c724d39350627ae Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Fri, 25 Sep 2026 14:09:20 +0200 Subject: [PATCH 5/8] WIP --- examples/python/scalar-field-mapping.py | 2 +- .../components/engine/GridSampler.h | 2 -- .../components/mapping/ScalarFieldMapping.cpp | 14 +++++++++----- 3 files changed, 10 insertions(+), 8 deletions(-) diff --git a/examples/python/scalar-field-mapping.py b/examples/python/scalar-field-mapping.py index b8f2730..b63272f 100644 --- a/examples/python/scalar-field-mapping.py +++ b/examples/python/scalar-field-mapping.py @@ -39,5 +39,5 @@ def createScene(root): model.state.showObject = True model.state.showObjectScale = 5.0 - model.Child.addObject("StopperLagrangianConstraint", name="constraint", min=0.0, max=1.0, index=1) + model.Child.addObject("StopperLagrangianConstraint", name="constraint", min=0.0, max=1.0, index=3) model.Child.addObject("GenericConstraintCorrection", name="correction", linearSolver=root.linearSolver.linkpath, ODESolver=root.odesolver.linkpath) diff --git a/src/SofaImplicitField/components/engine/GridSampler.h b/src/SofaImplicitField/components/engine/GridSampler.h index 1fcf5a1..ef1e4db 100644 --- a/src/SofaImplicitField/components/engine/GridSampler.h +++ b/src/SofaImplicitField/components/engine/GridSampler.h @@ -52,8 +52,6 @@ class GridSampler : public BaseComponent Data d_buffer; - void notifyLinkSet(BaseLink*, Base*) override; - protected: SingleLink l_field; diff --git a/src/SofaImplicitField/components/mapping/ScalarFieldMapping.cpp b/src/SofaImplicitField/components/mapping/ScalarFieldMapping.cpp index 3ffd076..3367157 100644 --- a/src/SofaImplicitField/components/mapping/ScalarFieldMapping.cpp +++ b/src/SofaImplicitField/components/mapping/ScalarFieldMapping.cpp @@ -110,15 +110,19 @@ void ScalarFieldMapping::applyJT( const ConstraintParams* mparams, InDataMatrixD if(!isComponentStateValid()) return; auto x = fromModel->readPositions(); + auto y = toModel->readPositions(); auto dx = sofa::helper::getWriteOnlyAccessor(dx_); auto dy = sofa::helper::getReadAccessor(dy_); auto field = l_field.get(); - for (Size i = 0; i < dx.size(); ++i) - { - const Vec3d grad = field->getGradient(x[i]); - dx[i] += grad * dy[i].x(); //< Because dy is in R so there is only one value - } + std::cout << "Constraint InDataMatrixDeriv:" << dx_ << std::endl; + std::cout << "Constraint OutDataMatrixDeriv:" << dy_ << std::endl; + + // for (Size i = 0; i < dx.size(); ++i) + // { + // const Vec3d grad = field->getGradient(x[i]); + // dx[i] += grad * dy[i].x(); //< Because dy is in R so there is only one value + // } } From 8186a2c8116d06210cf730402e645a07ec8291ed Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Fri, 25 Sep 2026 14:11:43 +0200 Subject: [PATCH 6/8] FIXUP --- src/SofaImplicitField/components/engine/GridSampler.cpp | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/SofaImplicitField/components/engine/GridSampler.cpp b/src/SofaImplicitField/components/engine/GridSampler.cpp index 1f4f16d..367e9b9 100644 --- a/src/SofaImplicitField/components/engine/GridSampler.cpp +++ b/src/SofaImplicitField/components/engine/GridSampler.cpp @@ -67,13 +67,6 @@ void GridSampler::init() d_componentState = core::objectmodel::ComponentState::Valid; } -void GridSampler::notifyLinkSet(BaseLink* link, Base* base) -{ - if(link==&l_field){ - std::cout << "THE FIELD HAS HCNAGD" << std::endl; - } -} - void GridSampler::updateInternalBuffer(const Vec3u& resolution, const Vec3d& min, const Vec3d& max) { auto buffer = sofa::helper::getWriteOnlyAccessor(d_buffer); From 425aa38673e8fa43d46b2db5d39dabbc802f1ed5 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Mon, 28 Sep 2026 12:26:58 +0200 Subject: [PATCH 7/8] Add IsoSurfaceSampler component Point point equally sampled on a given iso surface. --- CMakeLists.txt | 2 + .../components/engine/IsoSurfaceSampler.cpp | 164 ++++++++++++++++++ .../components/engine/IsoSurfaceSampler.h | 70 ++++++++ 3 files changed, 236 insertions(+) create mode 100644 src/SofaImplicitField/components/engine/IsoSurfaceSampler.cpp create mode 100644 src/SofaImplicitField/components/engine/IsoSurfaceSampler.h diff --git a/CMakeLists.txt b/CMakeLists.txt index b1dce2a..727026b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,6 +22,7 @@ set(HEADER_FILES ${SOFAIMPLICITFIELD_SRC}/components/engine/FieldToSurfaceMesh.h ${SOFAIMPLICITFIELD_SRC}/components/engine/GridSampler.h + ${SOFAIMPLICITFIELD_SRC}/components/engine/IsoSurfaceSampler.cpp ${SOFAIMPLICITFIELD_SRC}/components/geometry/BottleField.h ${SOFAIMPLICITFIELD_SRC}/components/geometry/DiscreteGridField.h ${SOFAIMPLICITFIELD_SRC}/components/geometry/SphericalField.h @@ -43,6 +44,7 @@ set(SOURCE_FILES ${SOFAIMPLICITFIELD_SRC}/components/engine/FieldToSurfaceMesh.cpp ${SOFAIMPLICITFIELD_SRC}/components/engine/GridSampler.cpp + ${SOFAIMPLICITFIELD_SRC}/components/engine/IsoSurfaceSampler.cpp ${SOFAIMPLICITFIELD_SRC}/components/geometry/BottleField.cpp ${SOFAIMPLICITFIELD_SRC}/components/geometry/ScalarField.cpp ${SOFAIMPLICITFIELD_SRC}/components/geometry/DiscreteGridField.cpp diff --git a/src/SofaImplicitField/components/engine/IsoSurfaceSampler.cpp b/src/SofaImplicitField/components/engine/IsoSurfaceSampler.cpp new file mode 100644 index 0000000..910ea86 --- /dev/null +++ b/src/SofaImplicitField/components/engine/IsoSurfaceSampler.cpp @@ -0,0 +1,164 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture, development version * +* (c) 2006-2025 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +****************************************************************************** +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include +#include + +#include +#include +using sofa::core::RegisterObject; +using sofa::core::visual::VisualParams; + +namespace sofaimplicitfield::component::engine +{ + +IsoSurfaceSampler::IsoSurfaceSampler() + : d_resolution(initData(&d_resolution, Vec1u{10}, "resolution", "Number of samples in each dimension")) + , d_min(initData(&d_min, {-1.0, -1.0, -1.0}, "min", "Minimum corner of the sampling grid")) + , d_max(initData(&d_max, {1.0, 1.0, 1.0}, "max", "Maximum corner of the sampling grid")) + , l_field(initLink("field", "The scalar field to sample")) +{ + addUpdateCallback("sample", {&d_resolution, &d_min, &d_max}, [this](const sofa::core::DataTracker&) + { + Vec1u resolution = d_resolution.getValue(); + Vec3d min = d_min.getValue(); + Vec3d max = d_max.getValue(); + + sampleField(); + return core::objectmodel::ComponentState::Valid; + }, {}); +} + +IsoSurfaceSampler::~IsoSurfaceSampler() +{ +} + +void IsoSurfaceSampler::init() +{ + if (!l_field.get()) + { + msg_error() << "Missing scalar field to sample"; + d_componentState = core::objectmodel::ComponentState::Invalid; + return; + } + + sampleField(); + d_componentState = core::objectmodel::ComponentState::Valid; +} + +void IsoSurfaceSampler::sampleField() +{ + auto field = l_field.get(); + if (!field) + { + msg_error() << "No scalar field linked"; + return; + } + + auto buffer = sofa::helper::getWriteAccessor(d_buffer); + auto resolution = buffer->resolution; + auto min = buffer->min; + auto max = buffer->max; + + if (resolution[0] <= 0 || resolution[1] <= 0 || resolution[2] <= 0) + { + msg_error() << "Resolution must be greater than 0 in all dimensions"; + return; + } + + if (min[0] >= max[0] || min[1] >= max[1] || min[2] >= max[2]) + { + msg_error() << "min must be less than max in all dimensions"; + return; + } + + std::cout << getPathName() << " sampling the grid " << std::endl; + auto data = buffer->data; + auto spacing = buffer->spacing; + + std::cout << " : " << buffer->data << std::endl; + + int rx = resolution.x(); + int rxy = resolution.y() * rx; + + std::vector positions; + std::vector results; + positions.reserve(rxy); + results.reserve(rxy); + + for (unsigned int z = 0; z < resolution[2]; z++) + { + positions.clear(); + results.clear(); + for (unsigned int y = 0; y < resolution[1]; y++) + { + for (unsigned int x = 0; x < resolution[0]; x++) + { + positions.emplace_back(min + Vec3d( + static_cast(x) * spacing[0], + static_cast(y) * spacing[1], + static_cast(z) * spacing[2] + )); + } + } + field->getValues(positions, results); + for (unsigned int y = 0; y < resolution[1]; y++) + { + for (unsigned int x = 0; x < resolution[0]; x++) + { + auto index=[rxy,rx](int x,int y, int z){ return z * rxy + y * rx + x; }; + auto rindex=[rx](int x,int y){ return y * rx + x; }; + + data[index(x,y,z)] = results[rindex(x,y)]; + } + } + } + std::cout << "SAMPLING DONE "<< std::endl; +} + +void IsoSurfaceSampler::computeBBox(const core::ExecParams* params, bool onlyVisible) +{ + if (onlyVisible) return; + + Vec3d min = d_min.getValue(); + Vec3d max = d_max.getValue(); + + sofa::core::objectmodel::BaseComponent::computeBBox(params, onlyVisible); + + f_bbox.setValue({min, max}); +} + +void IsoSurfaceSampler::draw(const VisualParams* vparams) +{ + if (!vparams || !vparams->displayFlags().getShowBehaviorModels()) + return; + + if (isComponentStateInvalid()) + return; +} + +void registerIsoSurfaceSampler(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData("Samples a ScalarField on a 3D grid and stores the result in a DiscreteGridField.") + .add< IsoSurfaceSampler >()); +} + +} diff --git a/src/SofaImplicitField/components/engine/IsoSurfaceSampler.h b/src/SofaImplicitField/components/engine/IsoSurfaceSampler.h new file mode 100644 index 0000000..a3aa4d5 --- /dev/null +++ b/src/SofaImplicitField/components/engine/IsoSurfaceSampler.h @@ -0,0 +1,70 @@ +/****************************************************************************** +* 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 + +//////////////////////////////////////////////////////////////////////////////////////////////////// +namespace sofaimplicitfield::component::engine +{ + +namespace{ + using namespace sofa; + using sofa::core::objectmodel::BaseComponent; + using sofa::core::visual::VisualParams; + using sofa::type::Vec3d; + using sofa::type::Vec1u; + using sofa::component::geometry::ScalarField; + using sofa::component::geometry::DiscreteGridField; +} + +class IsoSurfaceSampler : public BaseComponent +{ +public: + SOFA_CLASS(IsoSurfaceSampler, BaseComponent); + + void init() override; + void draw(const VisualParams* params) override; + + Data d_resolution; + Data d_min; + Data d_max; + + Data d_buffer; + +protected: + SingleLink l_field; + +protected: + IsoSurfaceSampler(); + virtual ~IsoSurfaceSampler(); + +private: + void computeBBox(const core::ExecParams* params, bool onlyVisible = false) override; + void updateInternalBuffer(const Vec1u& resolution, const Vec3d& min, const Vec3d& max); + void sampleField(); +}; + +} + From 02d3f818f755b2dfd9bc48c8c858200b7f298f83 Mon Sep 17 00:00:00 2001 From: Damien Marchal Date: Mon, 28 Sep 2026 12:27:53 +0200 Subject: [PATCH 8/8] Add test example. --- examples/python/scalarfield-forcefield.py | 41 +++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 examples/python/scalarfield-forcefield.py diff --git a/examples/python/scalarfield-forcefield.py b/examples/python/scalarfield-forcefield.py new file mode 100644 index 0000000..63d3b14 --- /dev/null +++ b/examples/python/scalarfield-forcefield.py @@ -0,0 +1,41 @@ +import Sofa +import Sofa.Core +from SofaImplicitField import ScalarField +import numpy + +def Model(): + particules = Sofa.Core.Node("Particles") + particules.addObject("MechanicalObject", name="state", template="Vec3", position=[[-1,0,0],[0,0,0],[1,0,0]]) + particules.addObject("UniformMass", name="mass", totalMass=1.0) + + child = particules.addChild("Child") + child.addObject("MechanicalObject", name="state", template="Vec1", position=[0,0,0]) + child.addObject("SphericalField", name="field", center=[0,0,0], radius=1.0) + child.addObject("ScalarFieldMapping", name="mapping", + input=particules.state.linkpath, + output=child.state.linkpath, field=child.field.linkpath) + + return particules + +def createScene(root): + """In this scene we create a ScalarField of spherical shape and a mechanical object of 3D particles. + The field and the particules are used as input of a ScalarFieldMapping, mapping the particules to a 1D child space. + This child space is then used in a constraints. + """ + root.addObject('RequiredPlugin', pluginName='Sofa.Component.AnimationLoop') # Needed to use components [FreeMotionAnimationLoop] + root.addObject('RequiredPlugin', pluginName='Sofa.Component.Constraint.Lagrangian.Model') # Needed to use components [StopperLagrangianConstraint] + root.addObject('RequiredPlugin', pluginName='Sofa.Component.Constraint.Lagrangian.Solver') # Needed to use components [BlockGaussSeidelConstraintSolver] + root.addObject('RequiredPlugin', pluginName='Sofa.Component.StateContainer') # Needed to use components [MechanicalObject] + root.addObject("RequiredPlugin", pluginName="SofaImplicitField") + root.addObject('RequiredPlugin', pluginName='Sofa.Component.Mass') # Needed to use components [UniformMass] + + root.addObject("DefaultAnimationLoop") + + root.addObject("EulerImplicitSolver", name="odesolver", rayleighStiffness=0.1, rayleighMass=0.1) + root.addObject("CGLinearSolver", name="linearSolver", template="CompressedRowSparseMatrixd") + + model = root.addChild(Model()) + model.state.showObject = True + model.state.showObjectScale = 5.0 + + model.Child.addObject("AllowedIntervalSpringForceField", name="forcefield")