Usage of volumetric meshes

Support
3

Hi Julian,

I thought I would share my experience that I have had lately with volumetric meshing and include the problems I am having

I tried using cleaver with limited success. I had a lot of trouble with the extension initially. After I spoke with Andras Lasso I was able to get the meshing working. The main thing I did to get cleaver to work with my data was to create a label map with a value of 75 on the inside and -75 on the outside. I used simple filters to do the shifting and scaling.

2 notes on using cleaver: 1) the mesh did not contain tetrahedra as far as I could tell. Instead it contained 4 triangles representing the faces instead of tetrahedra. 2) the mesh was not aligned in physical space with the volumes.

The method I am currently working with is using vtkDataSetTriangleFilter with something like the following at the python interpreter

>>> import vtk
>>> dstFilter = vtk.vtkDataSetTriangleFilter()
>>> croppedScan = getNode("croppedScan")
>>> dstFilter.SetInputData(croppedScan.GetImageData())
>>> dstFilter.Update()
>>> mesh = dstFilter.GetOutput()
>>> modelNode = slicer.vtkMRMLModelNode()
>>> modelNode.SetAndObserveMesh(mesh)
>>> slicer.mrmlScene.AddNode(modelNode)

This creates a solid rectangular prism mesh, where each voxel is split into several tets. The nodes are labeled with the intensity of the scan.I was able to cut out regions of the mesh based on the intensity at the nodes using vtkThreshold.

I am having trouble labeling the tetrahedra with intensity. I wrote the following function to label the tetrahedra with intensity. After this processing step the meshes show up as having zero cells in the Models module. Further they are not visualized in the 3D scene. The vtk files written from these models appear to have the cell data. Does anyone have any thought on how to go about labeling tets in Slicer that allows visualization? Is there some problem with my approach?

def sampleNodeScalarOntoCellScalar(mesh):
    """

    :type mesh: vtkUnstructuredGrid
    """
    idList = vtk.vtkIdList()
    cells = mesh.GetCells()
    nCell = cells.GetNumberOfCells()

    scanIntensity = vtk.vtkFloatArray()
    scanIntensity.SetNumberOfComponents(1)
    scanIntensity.SetNumberOfValues(nCell)

    for cellIndex in range(nCell):
        mesh.GetCellPoints(cellIndex, idList)
        scalarArray = mesh.GetPointData().GetArray(0)
        meanScalar = 0
        numPoints = idList.GetNumberOfIds()
        for pointIndex in range(numPoints):
            meanScalar = meanScalar + scalarArray.GetValue(idList.GetId(pointIndex))


        meanScalar = meanScalar / numPoints
        scanIntensity.SetValue(cellIndex, meanScalar)

    mesh.GetCellData().SetScalars(scanIntensity)
    return scanIntensity

Thanks for the help!

Michael

1 Like