Biomech GenIso Mesher

This mesher is used to generate a 3D surface mesh from a 3D segmented image.

The input is a 3D segmented binary image. It contains two values: voxels with a $0$ value are located outside of the object to mesh, voxels with a non-zero value are located inside.

The image must be .vtk, or event better .vti, which is more compact. If it is not the right format use Slicer3D, Paraview or my interface geniso/tools/gui4.py to change it.

First, the geniso tools must be loaded:

 import geniso.tools.generalTools as generalTools
 import geniso.tools.meshingTools as meshingTools
 import geniso.tools.imagingTools as imagingTools

Then the image must be loaded:

 image = generalTools.loadVtkImageXML(imageFileName)

Afterwards, the following set of commands is used to generate a surface mesh of the boundary of the segmented image:

 gen = GenerateMesh(generalTools.vtkImageDataToCpp(image)) 
 gen.execute() 
 smooth = SmoothMesh(gen.getSurf()) 
 smooth.smoothWithFctFromPoly(nbr,spx,spy,spz,nmin,interpol,type)

where the parameters are described in the table below.

In practice, GenerateMesh(image).execute() creates a triangulation of the boundary of the segmented object in the binary image given. Because of the spatial discretization of the image, the mesh which is generated is not smooth, since the steps of the image are kept. The class SmoothMesh() is used to eliminate these irregularities while keeping the geometry described in the image. It does so with the construction of an implicit function, described by the parameters nmin, interpol and type (for more details: http://hdl.handle.net/2268/136159).

Then, the vtkPolydata is retrieved with the command:

 poly = generalTools.genisoMeshToPolyData(smooth.getSurf())  

TetGen is used to get a tetrahedral volume mesh of the object from the triangular surface mesh of the boundary of the object:

 ugrid = meshingTools.callTetgen(poly)

A vtkUnstructuredGrid is retrieved. To save it on disk: generalTools.saveUgridXML('ugrid.vtu',ugrid). Paraview is used to view it. The class LoadMesh, implemented in the file, is used to include it in a Metafor test case:

 groupset.add(Group(1))
 mesh = LoadMesh(domain, ugrid, groupset(1) )
 mesh.execute()

This option is used to generate a mesh which contains several material areas: two different bones, a brain with ventricles separated from a tumor, a jawbone with a series of teeth…

In the multi-material case, the segment image representing the geometry must contain several labels: 0 for the background (outside of the object), and a series of non-zero positive values (typically 1, 2, 3…) for each material.

To use the GenIso mesher, the implicit function Fct must first be defined. This function is a continuous equivalency of the discrete segment image, and will allow the mesher to supply a smooth interface between the mesh separating each region.

To define this function, couples of labels of the image are defined, labels between which the mesher must create a triangulation. In addition, these couples are sorted to define the heterogeneous object as a set of closed surfaces on which open surfaces can be attached. In the image below, a tumor (open surface defined by the couple (4, 3)), is tied on the external boundary of the brain (closed surface defined by the union of the couples (3, 0) and (4, 0)).

Consequent, the syntax is as follow:

  fct = Fct()
  surf1 = fct.defineClosedSurface()
  surf1.mesh(3,0)
  surf1.mesh(4,0)
  surf2 = fct.defineOpenSurface()
  surf2.mesh(4,3)
  fct.setEpsilon(eps)
  fct.setNmin(nmin)
  fct.setType(type)
  fct.setInterpol(interpol)
  fct.setVtkImageFromPython(generalTools.vtkImageDataToCpp(image))
  fct.build()

Afterwards, the mesh can be generated with the commands:

  gen = GenerateMesh(fct)
  gen.getRegularGrid().setResolution(res)
  gen.execute()
  poly = generalTools.genisoMeshToPolyData(gen.getSurf()) // vtkPolyData

From this surface multi-regions mesh, TetGen is used to generate the volume mesh in each of these regions. To do so, TetGen must be given a point, seed situated in each of the regions to be meshed (this point can be found, for example, by loading the vtkPolydata in the graphical interface geniso/tools/gui4.py using the option point).

  seedbrain= [118.3,160.5,75]
  seedtumor = [118,196.9,71]
  regionSeeds = [seedbrain,seedtumor]
  regionVolmax = [volmax, volmax]
  ugrid = meshingTools.callTetgenMultipleRegions(poly, regionSeeds, regionVolmax)

Then, the mesh is loaded in Metafor:

  labels = [1,2]  // the two labels of the regions located in the ugrid (see visu in Paraview)
  grp = [groupset(1),groupset(2)] 
  
  mesh = tetgen.LoadMesh(domain, ugrid, grp)
  mesh.setFillMeshPoints(True)
  mesh.setGrpLabels(labels)
  mesh.execute()