I am working on a project to map material properties from CT scans to elements. I have been working on a couple ideas. The goal is to get the density of each tetrahedral element for a mesh of a bone that I got from segmenting in 3D Slicer. The goal is to get a good mix of accuracy and computation time.
The easiest one is to select the HU values at each node and average them, then divide by the volume of the tetrahedral to get the density. This has accuracy issues because a lot of voxels are not selected. I’ve coded and tested this already.
Next, I can use a lot of integration points to sample multiple HU values inside the tetrahedral in addition to the nodes of the tetrahedral. This is more computationally expensive but captures more voxels. I’ve coded and tested this already.
Lastly, segment the region of each tetrahedral over the CT scan and find the average HU value. I have not been able to code and automate the segmentation of this part by when providing the four nodes (help here; I couldn’t find relevant examples). This seems by far the slowest option because I would have to segment tens of thousands of tetrahedrals. There’s no way this is computationally efficient.
I use integration fields through integration points from trilinear interpolation. I am not familiar with the math but I do know at the integration points and just integration from the 4 points of the tetrahedral. There has been a study on this and is how Bonemat works. However, I have no clue how to implement this. The math seems hard to implement. I even looked at the python version of Bonemat and I can’t seem to figure it out. Paper: https://pubmed.ncbi.nlm.nih.gov/14644599/
If you guys have experience with this, let me know! I’d love to learn what’s best and what to implement.
You could create a scalar field with the element IDs and use the vtkProbeFilter with the CT grid (or subsampled if you want for efficiency) so you would get the id at each voxel. Then you can easily integrate.
Keep in mind though that if your elements are fairly large compared to the internal structure of the bone then simply averaging the Hu values won’t account for any differences in internal structure. You may want to include other features, like the variance or some frequency measure. But everything is an approximation anyway, so it would be good to have some experimental data to use for comparison of whatever technique you choose.
I was able to get the points from my mesh and the volume data; however, I don’t see a scripting example that I could use to extract HU from a tetrahedral region of interest. The closest examples I see use functions that support a cube or a sphere. Do you have some advice?
Here’s an example of using a vtkProbeFilter to extract a regular grid of displacements from a tetrahedral mesh. Extracting the cell IDs would be similar.
And here’s an example of getting the ct value at the centroid of tetrahedra to define Young’s modulus:
Here’s how you create cell scalars array:
So you would need to mix these together so you would get a labelmap of the cell IDs at each pixel and then you could map that info into your Young’s modulus map.
Let us know how it goes for you. These scripts are just experiments so far, but at somepoint we’ll want to have some useful high-level functionality for this kind of thing.
This is my current script. It’s able to take coordinates of nodes from a numpy array and get the voxel value at any coord point I input. This includes centroids and other integration points.
# import slicer
import slicer
import numpy as np
import csv
import time
nodes_coords = np.load(r"\node_Array.npy")
# Get your bone mask segment
volumeNode = slicer.util.getNode('Bone_Scan')
for row in nodes_coords:
# Get point coordinate in RAS
x = row[0] * -1
y = row[1] * -1
z = row[2]
point_Ras = [x,y,z]
# If volume node is transformed, apply that transform to get volume's RAS coordinates
transformRasToVolumeRas = vtk.vtkGeneralTransform()
slicer.vtkMRMLTransformNode.GetTransformBetweenNodes(None, volumeNode.GetParentTransformNode(), transformRasToVolumeRas)
point_VolumeRas = transformRasToVolumeRas.TransformPoint(point_Ras)
# Get voxel coordinates from physical coordinates
volumeRasToIjk = vtk.vtkMatrix4x4()
volumeNode.GetRASToIJKMatrix(volumeRasToIjk)
point_Ijk = [0, 0, 0, 1]
volumeRasToIjk.MultiplyPoint(np.append(point_VolumeRas,1.0), point_Ijk)
point_Ijk = [ int(round(c)) for c in point_Ijk[0:3] ]
voxelValue = voxels[point_Ijk[2], point_Ijk[1], point_Ijk[0]] # note that numpy array index order is kji (not ijk)
# print(row)
# print(point_Ras)
# print(point_Ijk)
# print(voxelValue)
I want to extract the voxel values within 4 points of a tetrahedral. Is there a method to apply a label map, segmentation, or ROI within the tetrahedral for me to extract the values?
I get this error when I try to run your first example.
typ = vtk_array.GetDataType()
AttributeError: 'NoneType' object has no attribute 'GetDataType'
[VTK] Input port 1 of algorithm vtkProbeFilter (000001F4B8128160) has 0 connections but is not optional.