Refactoring of DiscreteGridField and add ScalarFieldMapping

CMakeLists.txt

@@ -13,37 +13,43 @@ 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
     ${SOFAIMPLICITFIELD_SRC}/deprecated/ImplicitSurfaceContainer.h    # This is a backward compatibility file toward ScalarField
     ${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
     ${SOFAIMPLICITFIELD_SRC}/components/geometry/ScalarField.h
     ${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
     ${SOFAIMPLICITFIELD_SRC}/initSofaImplicitField.cpp
     ${SOFAIMPLICITFIELD_SRC}/MarchingCube.cpp
+    ${SOFAIMPLICITFIELD_SRC}/MHD.cpp
 
     ## This is a backward compatibility..
     ${SOFAIMPLICITFIELD_SRC}/deprecated/SphereSurface.cpp
     ${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
     ${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

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)
+

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=3)
+    model.Child.addObject("GenericConstraintCorrection", name="correction", linearSolver=root.linearSolver.linkpath, ODESolver=root.odesolver.linkpath)

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

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,15 +129,15 @@ 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;
 
         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;
@@ -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<ScalarField>([](sofa::core::objectmodel::Base* object) {
+        return py::cast(dynamic_cast<ScalarField*>(object));
+    });
 }
 
 }

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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: [email protected]                             *
+******************************************************************************/
+#include <SofaImplicitField/config.h>
+#include <SofaImplicitField/MHD.h>
+#include <sofa/type/Vec.h>
+#include <fstream>
+#include <cstring>
+#include <string>
+
+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;
+}
+
+
+}