I am trying to get iMin, iMax, jMin, jMax, kMin, kMax from some ROI node using Python script such that croppedVolume = dataVolume[kMin:kMax, jMin:jMax, iMin:iMax] and then I will perform some analysis on the cropped volume. What I know is that I can get the center and radius of ROI box in RAS coordinates by roiNode.GetXYZ() and roiNode.GetRadiusXYZ(). So I can get the RAS coordinates of center points of the 6 faces of the ROI box. Assuming the edges of ROI box are parallel to the RAS coordinate axes (which might not be true), I should be able to compute iMin, iMax, jMin, jMax, kMin, kMax with the help of dataVolumeNode.GetRAStoIJKMatrix(). For example, the converted ijk coordinates of the center points on left and right plane should have two same components and the differing component should give you iMin/iMax (or jMin/jMax or kMin/kMax). However, all three components I get are different. Then I saw the source code of VolumeClipWithROI here (starting from line 254), so I realized there might be another coordinate transform between ROI coordinates and volume RAS coordinates. But I get lost from there. So my questions are:
Is there a way to compute a slice object (in the form of iMin, iMax, jMin, jMax, kMin, kMax or other) such that croppedVolume = dataVolume[sliceObject]? BTW, I don’t care about the interpolation at the boundary. So if the resulting ijk are not integers, rounding is perfectly fine with me.
What if the ROI box is rotated such that the edges of the box are no longer parallel with the RAS coordinate axes? I suppose a binary mask might be needed to indicate which portion of the original volume gets cropped?
You can easily get any point’s position in any coordinate system by computing the transform between the coordinate systems (4x4 matrix) and multiply the point coordinates by the matrix.
The transformation chain: ROI -> RAS (a.k.a. World) -> Volume RAS -> Volume IJK. You can get ROI to Volume RAS directly using GetTransformBetweenNodes and append the volumes RAS to IJK transform. See some information and examples here:
Points to transform: Probably the best is to transform the 8 corner points of the ROI (you compute the 8 point coordinates in ROI coordinate system and transform one by one to Volume IJK).
What would you like to use the transformed ROI for? For images, it is often more efficient to crop the image with the ROI (to reduce the amount of data) - for example, using Crop volumes module.
If you need to define an arbitrarily shaped ROI, then you need to use image masks. You can create image masks using Segment editor module and export it to binary labelmap using Segmentations module.
GetTransformBetweenNodes seems to be what I need! I saw one of your earlier post on GetTransformBetweenNodes here, so my code should probably look like this:
transformMatrix = vtk.vtkMatrix4x4() # I think transforms between these coordinate systems should be linear?
a = getNode('myROI').GetParentTransformNode()
b = ?
slicer.vtkMRMLTransformNode().GetMatrixTransformBetweenNodes(a, b, transformMatrix)
ijkCoordinates = transformMatrix.MultiplyPoint(rasCoordinates of ROI points appended by 1)
But how should I define b? Since it represents the ijk coordinates to index the data matrix and it does not have a corresponding node in Slicer…
To your second question, the reason why I didn’t want to use Crop volumes module is that I will do some calculation on the cropped volume, change its value, and eventually put it back to the original volume (so I need to get the ijk coordinates to know where to put it back)
a = ROI node
b = volume node
Once you have the matrix between a and b, multiply it with Volume RAS to IJK matrix (you can get this matrix from the volume node) to compute ROI to IJK. Pay attention to correct ordering (multiplication is not commutative for matrices) and if you need to invert matrices.
Cropping does not change physical position of voxels, so you can always copy the contents back to the original volume.
@lassoan I updated my code and now I get volume RAS coordinates of the center point in the format of (a1,b,c),(a2,b,c),(a,b1,c),(a,b2,c),(a,b,c1),(a,b,c2), which is as expected. However, when I applied RAS to IJK transformation, none of the three components are the same anymore. So I guess the volume RAS coordinates are not aligned with data matrix ijk axes and so I cannot use croppedVolume = dataVolume[kMin:kMax, jMin:jMax, iMin:iMax].
From your description, cropping seems to perfectly meet my needs. However, I couldn’t find any instructions on how to run the cropping module in python except this ClipVolumeWithRoi script. I am wondering if there is a way to do things like this:
croppedVolumeNode = crop(volumeNode, roiNode)
Perform some calculation on cropped volume...
volumeNode.update(croppedVolumeNode) # put it back to the original volume
Thank you. But may I ask how to put the cropped volume back? I see that the size of the cropped volume changes but the spatial location doesn’t. So is there a function I can use to replace the part of the original volume with the modified cropped volume? Thanks!
There are many ways to do this. If you have a rectangular ROI then you can use SimpleITK’s PasteImageFilter (it is available in Simple Filters module with a GUI, so you can find out how to use it). If you have an image mask, then create two images where parts to be replaced are set to a very small value then combine them using maximum operation (using numpy, ITK, or VTK filters, or Slicer modules).
I am sorry to bother you again on this but I tried the Simple Filters module and I am confused about DestinationIndex. I read the documentation on the PasteImageFilter and thought it should be one of the 8 corner coordinates of the ROI box expressed in volume RAS coordinates. I tried all 8 of them but none of them works (not even close). Could you please elaborate on how to compute DestinationIndex from the roiNode and volumeNode? BTW, I am using rectangular ROI hand dragged in 3D slicer.
Thanks!
This filter operates in index space. I’d start by leaving the SourceIndex, and Destination index as (0,0,0), then just specify the Source size. This will result in the filter just copying one corner from the source image to the destination image. Make sure you create a new output volume and select the proper input volumes.
Hello @blowekamp, I figured out the DestinationIndex issue, it was because when I created a new volume node in Python, I forgot to set the IJKToRASMatrix of it to be the same as that of the original volume. But now I have another problem here: after I put back the cropped volume back to the original volume, it seems that there is a small offset (looks like 1 pixel?) as below:
Using the crop volume module to crop the volume. I have also tried using the following code with the same result: slicer.modules.cropvolume.logic().CropVoxelBased(getNode('AnnotationROI_30'), getNode('MRHead'), croppedVolumeNode)
where croppedVolumeNode is some scalar volume node I created using scripts.
In the simple filters module, set the first input volume as MRHead, the second as MRHead cropped, source size is the same as the size of the cropped volume. Source and destination index are kept as (0,0,0). I also chose to create a new volume as the output.
I am wondering if I have done any steps wrong that produce this slight offset. And also is there a way to programmatically invoke the PasteImageFilter() in Python? I read the documentation page and found out the following code is needed:
myFilter = PasteImageFilter()
myFilter.SetDebug(False)
myFilter.SetDestinationIndex((0, 0, 0))
myFilter.SetNumberOfThreads(8)
myFilter.SetSourceIndex((0, 0, 0))
myFilter.SetSourceSize((194, 80, 129))
myFilter.SetSourceImage(??) # Should I put node or array here?
myFilter.SetDestinationImage(??)
But obviously this gives an error saying that PasteImageFilter() is not found. I also tried to import itk by using import itk but returns a module not found error. So could you please help me complete this snippet? Thanks!
So the SimpleFilters module utilizes SimpleITK for python. The interface is very similar to ITK’s but is designed to have more of a native Python feel, and compatible with may Python intrinsic types. You can find the SimpleITK documentation for the PasteImageFilter here. SimpleITK also has a procedural interface. Yesterday I did notice that the image inputs were not named. The first argument is the “DestinationImage”, while the second is the “SourceImage”. I have apull request that will improve this.
Your simpleITK code would look something like this:
Regarding your off by one error. I would closing inspect the image’s meta data, and the coordinations of the ROI. There may be slight differences in the physical locations of the image’s pixels that are causing problems, or look for rounding issues.
