How to get the IJK coordinate of a fiducial node in 3d slicer development

[zhutouzjq]

thank the help provided by jamesobulter and cpinter, now i can successfully choose the fiducial in myself extension.

but i fail to get the ijk coordinate of this fudicial node in my extension. i can only get the ras coordinate through the python code in class called ‘Logic‘.

def process(self,

            inputVolume: vtkMRMLScalarVolumeNode,

            imageThreshold: float,

            inputlowerpoint, #this is the name of qMRMLNodeCombobox by which i choose the fiducial node in 3D slicer

            inputupperpoint,

            showResult: bool = True) -> None:

    ras = vtk.vtkVector3d(0,0,0)

    inputlowerpoint.GetNthControlPointPositionWorld(0, ras)

    print('lowerpoint : IJK =’,ras)

after running the extension, the result in pyhton console like this:

however, what i need is the IJK coordinate of this fudicial node. So, how can i get the IJK coordinate of a fiducial node and print in the python console?

[muratmaga]

AI is great for this type of things. I asked your question to Gemini:

“how can I get the IJK coordinates of Fiducial point that I know its world coordinates (RAS) using Slicer’s python API? Assume the coordinate is in the first control point of PointList F”

and I got

import slicer
import vtk
import numpy as np

# 1. Define the RAS coordinates (World Coordinates)
ras_point = np.array([-38.0, 24.0, -617.0, 1.0]) # Homogeneous coordinates (R, A, S, 1)

# 2. Get the reference volume node
# Replace 'YourVolumeNodeName' with the actual name or ID of your loaded volume node.
# For example, if you have the built-in MRHead:
volumeNode = slicer.util.getNode('MRHead')

if volumeNode is None:
    print("Error: Volume node not found. Please load a volume into the scene.")
else:
    # 3. Get the WorldToIJK transform matrix from the volume node
    worldToIjkMatrix = vtk.vtkMatrix4x4()
    volumeNode.GetWorldToIJKMatrix(worldToIjkMatrix)
    
    # You can convert the vtkMatrix4x4 to a numpy array for easier calculation
    worldToIjkNp = slicer.util.vtkMatrixToNumpy(worldToIjkMatrix)
    
    # 4. Calculate the IJK coordinates
    # Matrix multiplication: [I, J, K, 1] = WorldToIJK_Matrix @ [R, A, S, 1]
    ijk_point_homog = worldToIjkNp @ ras_point
    
    # The IJK coordinates are the first three components (the last is 1.0)
    # The coordinates are typically rounded to the nearest integer for voxel indices.
    ijk_point = np.round(ijk_point_homog[:3]).astype(int)
    
    print(f"RAS Coordinates (R, A, S): {ras_point[:3]}")
    print(f"IJK Coordinates (I, J, K): {ijk_point}")
    
    # Optional: Display the IJK coordinates for a specific fiducial point
    # in an existing markups node (if needed for context)
    # fiducialNode = slicer.util.getNode('F') # Replace 'F' with your markups node name
    # if fiducialNode:
    #     fiducialNode.SetNthControlPointPositionWorld(0, ras_point[0], ras_point[1], ras_point[2])
    #     # After setting the position, you can use the node's internal function to get IJK:
    #     # ijk = [0, 0, 0]
    #     # fiducialNode.GetNthControlPointPositionIJK(0, ijk)
    #     # print(f"IJK from Markups node: {ijk}")

which seems to work.

[zhutouzjq]

thank the discussion in “RAS to IJK conversion issues“. the reply by muratmaga help me a lot.

in the script repository of 3D slicer, they suggest use the follow code to transform the ras to ijk:

# Inputs
volumeNode = getNode("MRHead")
pointListNode = getNode("F")
markupsIndex = 0

# Get point coordinate in RAS
point_Ras = [0, 0, 0]
pointListNode.GetNthControlPointPositionWorld(markupsIndex, point_Ras)

# 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] ]

# Print output
print(point_Ijk)

but what i didn’t know is that in my code, the second line “volumeNode = getNode(“MRHead”)” should be change into “volumeNode = slicer.util.getNode(“5: CTA 0.5 CE”)“.

“MRHead“ or “5: CTA 0.5 CE“ is the name of volume which could be found here.

and use the python code like this in Logic class, i can print the IJK coordinate in console now.

def process(self,

            inputVolume: vtkMRMLScalarVolumeNode,

            imageThreshold: float,

            inputlowerpoint,

            inputupperpoint,

            showResult: bool = True) -> None:

    """

    Run the processing algorithm.

    Can be used without GUI widget.

    :param inputVolume: volume to be thresholded

    :param outputVolume: thresholding result

    :param imageThreshold: values above/below this threshold will be set to 0

    :param invert: if True then values above the threshold will be set to 0, otherwise values below are set to 0

    :param showResult: show output volume in slice viewers

    """



    if not inputVolume:

        raise ValueError("Input volume is invalid")



    import time



    startTime = time.time()

    logging.info("Processing started")



    volumeNode = slicer.util.getNode("5: CTA 0.5 CE")

    markupsIndex = 0

    point_Ras = \[0, 0, 0\]

    inputlowerpoint.GetNthControlPointPositionWorld(markupsIndex, point_Ras)

    

    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\] \]

    print('lowerpoint : RAS =',point_Ras)

    print('lowerpoint : IJK =',point_Ijk)



    stopTime = time.time()

    logging.info(f"Processing completed in {stopTime-startTime:.2f} seconds")

[zhutouzjq]

thank you for your reply. with your discussion under the issue “RAS to IJK conversion issues“, i have solve this problem. your are very professional in developping 3d slicer extension