Coordinates of pixels in a mesh file

Support
1

Hello everyone,

I’m trying to get the coordinates of each voxel in my model through the following code:

volumeNode=slicer.util.getNode(“myfile.nrrd”)
ijkToRas = vtk.vtkMatrix4x4()
imageData=volumeNode.GetImageData()
extent = imageData.GetExtent()
for k in range(extent[4], extent[5]+1):
for j in range(extent[2], extent[3]+1):
for i in range(extent[0], extent[1]+1):
position_Ijk=[i, j, k, 1]
position_Ras=ijkToRas.MultiplyPoint(position_Ijk)
print(position_Ras)

The position_Ras are in the range of (-61.955299377441435, -57.42440032958985, 2.3572499752044687, 1.0) to (-13.861799377441436, -19.522900329589852, 21.785749975204467, 1.0) with Spacing: [0.3185, 0.3185, 0.3185] and Dimensions: (151, 119, 61).

So far, everything looks fine. However, when I look at my mesh generated using exported stl file in Gmsh, the node coordinates are in the range of (15.73733139, 21.27777481, 2.123630524) to (61.77390289, 56.87292862, 21.68377686).

Does anyone know how I can correlate these two coordinates? I need to specify which voxels are inside of an element in my mesh.

Thank you for your help,
Zahra

2

Hello @ZSoltani, the ijkToRas matrix does not seem to be set properly (it is only set to default values of the vtkMatri4x4 at creation). vtkMRMLVolumeNode has the GetIJKToRASMatrix() that you can use to populate the ijkToRas matrix.

If I understand correctly, you are comparing the extent of an image with the extent of a 3D surface model derived from the image. Most likely, the model does not occupy the whole volume of the image, but a subvolume of it, so your values are probably not going to match.

3

Also note that coordinates in all files (image files, mesh files, etc) are stored in LPS coordinate system by default. In Slicer, coordinates are in RAS coordinate system.

4

Thank you @RafaelPalomar and @lassoan for your replies.

Using volumeArray = slicer.util.array(“my-stl-file.stl”) I can get the same node coordinates as what I get from Gmsh. The only difference is that x and y coordinates have opposite signs in 3Dslicer and Gmsh and z coordinates are exactly the same. Is it a general rule? And if not, how I can transfer between two coordinate systems for a general stl file? Is it like they have different origins?

My second question, is about finding the coordinates of voxels inside the stl model. Sorry I’m new to 3Dslicer, and not sure how I can apply your comment correctly. Could you please explain with a short script how to do that?

Thank you,
Zahra

5

Here’s a description of LPS/RAS: Coordinate systems - Slicer Wiki

Yes, it’s just a matter of knowing the convention of the software you are using and flipping the signs on the first two coordinates to convert.

For your second question these scripts should help: https://slicer.readthedocs.io/en/latest/developer_guide/script_repository.html#get-volume-voxel-coordinates-from-markup-control-point-ras-coordinates

6

Thank you Steve for your help.