I have a list of RAS points I would like to assign to a pre-existing segment, either directly or via a label map. I looked at what’s in the nightly script repository, and it doesn’t appear to cover this situation.
I like the design of the segmentation nodes. I noticed that the label map associated with a segment has a minimal size. If there is only one labelled voxel, the array associated with the segment’s label map will only have one element.
Given this, I don’t know how to add a new point to a segment. I assume the array needs to be resized, and a new transform assigned. Is there a method that accomplishes this? How do the paint tools accomplish this?
Thanks.
If you implement this in a segment editor effect then you can use modifySelectedSegmentByLabelmap method (most effects use this method to modify segments).
If you need to modify a segment from outside of a segment editor effect, then probably the easiest is to add a new segment (that has an image data with the additional voxels painted) and add it to the existing segment by using Logical effect. Or, just use vtkImagePadFilter to change extents of your image data (no need to change origin or apply any transform, as only the image extents will be changed; having negative extents are completely fine).
1 Like
In the second case @lassoan mentions, instead of the vtkImagePadFilter you can use the convenience function vtkOrientedImageDataResample::PadImageToContainImage if imakes it easier.
2 Likes
Thank you @lassoan and @cpinter for your help. I’ve made some progress, but could still use some assistance.
Ultimately, what I’m doing is creating a new segmentation based on another segmentation, so the filtered one is fully contained in the range defined by the original. The original one is created using editor effects, while the filtered one is intended to appear automatically over it. Obtaining the points to assign to the filtered segmentation as a list of RAS points may not have been the best approach, but it’s the path I started down.
The first question I have is about the result: I get a grouping of points with individual boundaries drawn around them. I don’t know how to turn the result into a segmentation with one overall boundary. An image is attached.
The second question I have is about the efficiency of my approach. Is there a better way to accomplish what this method does?
Thanks again.
def NewAssignPointsToSegmentation(self, listOfPoints, originSegNodeID, originSegNumber, destinationSegNodeID, destinationSegNumber):
"""Assign a value at the given points in the segment with the given ID in the destination segmentation node."""
if(len(listOfPoints)==0):
#Clear the destination segment
return self.ClearSegmentation(destinationSegNodeID, destinationSegNumber)
#Get origin node
originSegNode = slicer.util.getNode(originSegNodeID) #type is vtkMRMLSegmentationNode
#... error checks
#Access the origin segment number specified in the parameter list
originSegmentation = originSegNode.GetSegmentation() #type is vtkSegmentation
originSegmentID = originSegmentation.GetNthSegmentID(originSegNumber)
#Get destination node
destinationSegNode = slicer.util.getNode(destinationSegNodeID) #type is vtkMRMLSegmentationNode
#... error checks
#Access the origin segment number specified in the parameter list
destinationSegmentation = destinationSegNode.GetSegmentation() #type is vtkSegmentation
destinationSegmentID = destinationSegmentation.GetNthSegmentID(destinationSegNumber)
#Segmentations module logic class (for static methods access)
segmentationsLogic = slicer.modules.segmentations.logic()
#Get the oriented image data for the originSegNode
originLabelMapVolumeNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLLabelMapVolumeNode')
segmentationsLogic.ExportAllSegmentsToLabelmapNode(originSegNode, originLabelMapVolumeNode)
originOrientedImage = segmentationsLogic.CreateOrientedImageDataFromVolumeNode(originLabelMapVolumeNode)
#Create the label node for the destination for the duration of this method
destinationLabelMapVolumeNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLLabelMapVolumeNode')
segmentationsLogic.ExportAllSegmentsToLabelmapNode(destinationSegNode, destinationLabelMapVolumeNode)
destinationOrientedImage = segmentationsLogic.CreateOrientedImageDataFromVolumeNode(destinationLabelMapVolumeNode)
#Increase the size of the destination image
vtkSegmentationCore.vtkOrientedImageDataResample.PadImageToContainImage(destinationOrientedImage, originOrientedImage, destinationOrientedImage)
#Fill destination image with zeros
vtkSegmentationCore.vtkOrientedImageDataResample.FillImage(destinationOrientedImage, 0)
#Get transform matrices
IJKToRASMatrix = vtk.vtkMatrix4x4()
originLabelMapVolumeNode.GetIJKToRASMatrix(IJKToRASMatrix)
RASToIJKMatrix = vtk.vtkMatrix4x4()
RASToIJKMatrix.DeepCopy(IJKToRASMatrix)
RASToIJKMatrix.Invert()
#Set new transform matrix for the destination oriented image
destinationOrientedImage.SetImageToWorldMatrix(IJKToRASMatrix)
valueToAssign = 1.0
#Add points to the destination oriented image
for labelPoint in listOfPoints:
rasPos = list(labelPoint[:3])
rasPos.append(1)
ijkPos = RASToIJKMatrix.MultiplyPoint(rasPos)
#Get IJK indices
i = int(round(ijkPos[0],0))
j = int(round(ijkPos[1],0))
k = int(round(ijkPos[2],0))
destinationOrientedImage.SetScalarComponentFromDouble(i,j,k,0,valueToAssign)
#Assign back to the segment given by destinationSegmentID
modeReplace = slicer.modules.segmentations.logic().MODE_REPLACE
destinationExtent = originLabelMapVolumeNode.GetImageData().GetExtent()
slicer.modules.segmentations.logic().SetBinaryLabelmapToSegment(destinationOrientedImage, destinationSegNode, destinationSegmentID, modeReplace, destinationExtent)
destinationSegNode.Modified()
#Clean up - delete the label map nodes
if(originLabelMapVolumeNode is not None):
slicer.mrmlScene.RemoveNode(originLabelMapVolumeNode)
if(destinationLabelMapVolumeNode is not None):
slicer.mrmlScene.RemoveNode(destinationLabelMapVolumeNode)
return True
#end NewAssignPointsToSegmentation
Do you know what may cause the boundary pixellation in the image I added?
Thanks.
Can you take a screenshot of the Segmentations module GUI to see what visualization settings you have there?
Ok. Here is the whole screen, with the Segmentations module at right. Do you need to see the Segment Editor GUI as well?
Display section in the Segmentations module widget is disabled. It either means that you don’t have a display node for your segmentation or there is something wrong with this widget. What do you see if you switch to Segmentations module? Is Display section disabled there, too? How does the segmentation looks in 3D and in other slice views?
You’re right, if I switch to the Segmentations module, the display section is enabled.
I don’t see the segmentations on 3D, and it is unable to convert any of them to closed surfaces.
It’s quite possible I’m doing something wrong with the display nodes. I am not sure how to handle those. I’ll post my segmentation configuration method.
def ConfigureSegmentations(self):
"""Configure the Segmentation nodes used in this module."""
SSFPSegmentationColor = self.ConfigureSegmentationColor('muscle')
UTETissueSegmentationColor = self.ConfigureSegmentationColor('tissue')
UTEFasciaOnlySegmentationColor = self.ConfigureSegmentationColor('connective tissue')
#Create the segmentation nodes
self.SSFPSegmentationNode = slicer.vtkMRMLSegmentationNode()
self.UTETissueSegmentationNode = slicer.vtkMRMLSegmentationNode()
self.UTEFasciaOnlySegmentationNode = slicer.vtkMRMLSegmentationNode()
slicer.mrmlScene.AddNode(self.SSFPSegmentationNode)
slicer.mrmlScene.AddNode(self.UTETissueSegmentationNode)
slicer.mrmlScene.AddNode(self.UTEFasciaOnlySegmentationNode)
self.SSFPSegmentationNode.CreateDefaultDisplayNodes()
self.UTETissueSegmentationNode.CreateDefaultDisplayNodes()
self.UTEFasciaOnlySegmentationNode.CreateDefaultDisplayNodes()
SSFPVolumeNode = self.GetSSFPVolumeNodeReference()
UTEVolumeNode = self.GetUTEVolumeNodeReference()
self.SSFPSegmentationNode.SetReferenceImageGeometryParameterFromVolumeNode(SSFPVolumeNode)
self.UTETissueSegmentationNode.SetReferenceImageGeometryParameterFromVolumeNode(UTEVolumeNode)
self.UTEFasciaOnlySegmentationNode.SetReferenceImageGeometryParameterFromVolumeNode(UTEVolumeNode)
SSFPSegmentation = self.SSFPSegmentationNode.GetSegmentation()
UTETissueSegmentation = self.UTETissueSegmentationNode.GetSegmentation()
UTEFasciaOnlySegmentation = self.UTEFasciaOnlySegmentationNode.GetSegmentation()
self.SSFPSegID = SSFPSegmentation.GenerateUniqueSegmentID('SSFPSeg')
self.UTETissueSegID = UTETissueSegmentation.GenerateUniqueSegmentID('UTETissueSeg')
self.UTEFasciaOnlySegID = UTEFasciaOnlySegmentation.GenerateUniqueSegmentID('UTEFasciaSeg')
verifySSFP = SSFPSegmentation.AddEmptySegment(self.SSFPSegID, 'sartorial', SSFPSegmentationColor)
verifyUTETissue = UTETissueSegmentation.AddEmptySegment(self.UTETissueSegID, 'tissuelabel', UTETissueSegmentationColor)
verifyUTEFasciaOnly = UTEFasciaOnlySegmentation.AddEmptySegment(self.UTEFasciaOnlySegID, 'fascia', UTEFasciaOnlySegmentationColor)
#end ConfigureSegmentationsWhat do you mean? How do you try and what happens instead? Is there an error message?
In the Representations box in the Segmentations module, if I click on Create, and then Convert in the dialog box that appears, I get a message box saying “Failed to convert Segmentation to Closed surface!”
Sorry, Create, then Update, then Convert.
I just noticed something interesting. Not all of the pixels are isolated. There are some pairs. In this slice, I can see two pairs.
You only need to click Create. Update is used for specifying conversion settings, such as smoothing. Do you do something like this?
The log file would help a lot to identify why conversion failed.
Also your data seems wrong overall. If you could send a scene (mrb) what would clear things up.
Ok. Am I able to attach it here, or should I email it to you?
The code used to create the fascia segmentation is in an earlier post. There are many things that could be wrong with how I’m creating it. If you can see a better way to do this, I’d be happy to fix it.
I got the scene via email, thanks! There are two problems:
- The scalar range of binary labelmap in the segment ‘fascia’ is not (0,1) as expected, but (0,254). I did a threshold like this to fix it (I was very sloppy with variable names etc, but you get the gist):
s=getNode('Segmentation_2')
im=s.GetBinaryLabelmapRepresentation('UTEFasciaSeg')
im.GetScalarRange() # Will be (0.0, 254.0)
t=vtk.vtkImageThreshold()
t.SetInputData(im)
t.ThresholdByUpper(1)
t.SetInValue(1)
t.Update()
im2=t.GetOutput()
im2.GetScalarRange() # Will be (0.0, 1.0)
im.DeepCopy(im2)
- The 3D viewer was turned off in the display node of this segmentation. I went to Display/Advanced/Views, and checked “View1” (only “Yellow” was checked)
I also disabled smoothing so that no data is lost from this sparse one-slice labelmap: when you click Update for Closed surface you see the list of conversion parameters, and by setting Smoothing to 0, it is disabled. You can do this programmatically by
s.GetSegmentation().SetConversionParameter('Smoothing factor', '0.0')
1 Like
If you import the labelmap to the segmentation and not just set it directly (I’m not sure how you do it), then it will be thresholded automatically:
vtkSlicerSegmentationsModuleLogic::ImportLabelmapToSegmentationNode
1 Like
Great, thanks! I’ll look in to this, and let you know how it goes.