Numpy array from slice, then get back to RAS coordinate system

Development
1

Hello,

Similar to this post (link) I am trying to:

  1. Visualize in the view node a slice that is perpendicular to a curve defined by two control points
  2. Generate a numpy array from slice (similar to the straighten volume example script)
  3. Process the numpy array to identify a single pixel coordinate (ijk) of interest
  4. Plot the ijk coordinate of the point of interest, in this case the middle of the array, on 3D slicer as marker (as in the post above). This new marker should align with the original marked used to get the slice.

However, my calculated coordinates are off center (see attached image). The point with the label “MarkupsCurve” is a manually placed center point used to get the slice and the point without a label is the point placed by the python script.

image
image1920×1059 110 KB

point_list = list(np.arange(numberOfPoints))
for pointIndex in point_list:
    print(pointIndex)
    curvePointToWorldMatrix = vtk.vtkMatrix4x4()
    curveNode.GetCurvePointToWorldTransformAtPointIndex(pointIndex, curvePointToWorldMatrix)
    curvePointToWorldTransform.SetMatrix(curvePointToWorldMatrix)
    sliceToWorldTransform.Update()
    sliceNode.GetSliceToRAS().DeepCopy(sliceToWorldTransform.GetMatrix())
    sliceNode.UpdateMatrices()
    slicer.app.processEvents()

    
    slicer.util.setSliceViewerLayers(background = volume)
    slicer.app.processEvents()
    slicer.util.forceRenderAllViews()
    # Get RAW slice and add it to the straightened volume
    tempSlice_RAW = vtk.vtkImageData()
    tempSlice_RAW.DeepCopy(reslice.GetOutput())
    append_RAW.AddInputData(tempSlice_RAW)
    
#     Get slice as volume to output into numpy array for processing by external network
    volumeNodeResliced_RAW.SetIJKToRASMatrix(sliceNode.GetXYToRAS())
    volumeNodeResliced_RAW.SetAndObserveImageData(tempSlice_RAW)
    volumeNodeResliced_RAW.CreateDefaultDisplayNodes()
    volumeNodeResliced_RAW.CreateDefaultStorageNode()
    voxels_RAW = slicer.util.arrayFromVolume(volumeNodeResliced_RAW)
    slices_RAW.append(voxels_RAW)
    
    # find RAS of middle pixel position in array - dummy example 
    print(voxels_RAW.shape)
    point_ijk = [voxels_RAW.shape[1]/2, voxels_RAW.shape[2]/2, 0]
    print(point_ijk)
    #Get physical coordinates from voxel coordinates
    volumeIjkToRas = vtk.vtkMatrix4x4()
    volumeNodeResliced_RAW.GetIJKToRASMatrix(volumeIjkToRas)
    point_VolumeRas = [0, 0, 0, 1]
    volumeIjkToRas.MultiplyPoint(np.append(point_ijk,1.0), point_VolumeRas)
    #If volume node is transformed, apply that transform to get volume’s RAS coordinates
    transformVolumeRasToRas = vtk.vtkGeneralTransform()
    slicer.vtkMRMLTransformNode.GetTransformBetweenNodes(volumeNodeResliced_MASK.GetParentTransformNode(), None, transformVolumeRasToRas)
    print(point_VolumeRas)
    point_Ras = transformVolumeRasToRas.TransformPoint(point_VolumeRas[0:3])
    print(point_Ras)
    curveNode.AddControlPoint(vtk.vtkVector3d(point_Ras))
    break
2

I would recommend to use the Curved Planar Reformat module to get the slices of the curve (or straight line). You can then process each slice of the straightened volume, get the physical (RAS) coordinates from voxel (IJK) coordinates, add the markup point using the RAS coordinates in each slice. You can apply the inverse of the straightening transform to move the points to the original image’s space.