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")
