|
| 1 | +#STEP 1 |
| 2 | +ID = "myskullid"##################type the name you want for the segmentation & model here |
| 3 | +mythresh = 2# input the threshold thickness value you want to dissplay on the colormap label |
| 4 | +#NOTE - this scripts interacts with items (volumes, segmentations, models) based on the name 'getNode('xxxx'), |
| 5 | +#so if you modify names, you will also need to modify the corresponding part of the script |
| 6 | + |
| 7 | +######################################################SETUP ENVIRONMENT########################################### |
| 8 | +#This section creates objects to access volumes, segmentations, and segmentation effects |
| 9 | + |
| 10 | +#set up the volume - name the volume as an object |
| 11 | +masterVolumeNode = slicer.mrmlScene.GetFirstNodeByClass("vtkMRMLScalarVolumeNode") |
| 12 | +volumeScalarRange = masterVolumeNode.GetImageData().GetScalarRange()#value used in threshold computation for segmentation |
| 13 | + |
| 14 | +# Create SEGMENTATION NODE AND LINK TO VOLUME |
| 15 | +segmentationNode = slicer.vtkMRMLSegmentationNode()#name segmentation node |
| 16 | +slicer.mrmlScene.AddNode(segmentationNode)#add the node to the scene |
| 17 | +segmentationNode.CreateDefaultDisplayNodes() # only needed for display |
| 18 | +segmentationNode.SetReferenceImageGeometryParameterFromVolumeNode(masterVolumeNode)#link the segmentation to the volume |
| 19 | + |
| 20 | +# Create segment editor to get access to segmentation effects |
| 21 | +segmentEditorNode = slicer.vtkMRMLSegmentEditorNode() |
| 22 | +slicer.mrmlScene.AddNode(segmentEditorNode)#add segment editor node to scene |
| 23 | +segmentEditorWidget = slicer.qMRMLSegmentEditorWidget() |
| 24 | +segmentEditorWidget.setMRMLScene(slicer.mrmlScene)#connect widget to scene |
| 25 | +segmentEditorWidget.setMRMLSegmentEditorNode(segmentEditorNode)#connect segment editor to editor widget |
| 26 | +segmentEditorWidget.setSegmentationNode(segmentationNode)#connect segmentation node |
| 27 | +segmentEditorWidget.setMasterVolumeNode(masterVolumeNode)#connect master node |
| 28 | + |
| 29 | +#calculate the threshold value for maximum entry algorithm segmentation |
| 30 | +import vtkITK |
| 31 | +ME_thresholdCalculator = vtkITK.vtkITKImageThresholdCalculator() |
| 32 | +ME_thresholdCalculator.SetInputData(masterVolumeNode.GetImageData()) |
| 33 | +ME_thresholdCalculator.SetMethodToMaximumEntropy() |
| 34 | +ME_thresholdCalculator.Update() |
| 35 | +Maxentval = ME_thresholdCalculator.GetThreshold() |
| 36 | + |
| 37 | +#STEP 2 |
| 38 | +####################################################Create skull segmentation###################################################### |
| 39 | +#create segmentation of the bone using the threshold effect and 'maximum entropy' algorithm |
| 40 | +skull = segmentationNode.GetSegmentation().AddEmptySegment(ID) |
| 41 | +segmentEditorNode.SetOverwriteMode(slicer.vtkMRMLSegmentEditorNode.OverwriteNone) |
| 42 | +segmentEditorNode.SetSelectedSegmentID(ID) |
| 43 | +segmentEditorNode.SetMaskSegmentID(ID) |
| 44 | +segmentEditorWidget.setActiveEffectByName("Threshold") |
| 45 | + |
| 46 | +#apply threshold effect |
| 47 | +effect = segmentEditorWidget.activeEffect() |
| 48 | +effect.setParameter("MinimumThreshold", str(Maxentval)) |
| 49 | +effect.setParameter("MaximumThreshold",str(volumeScalarRange[1])) |
| 50 | +effect.self().onApply()#apply separate |
| 51 | +segmentationNode.GetSegmentation().GetSegment(ID).SetColor(0,0,1) #color of segmentation is set to blue here |
| 52 | + |
| 53 | +#apply 'keep largest' island segmentation effect |
| 54 | +segmentEditorWidget.setActiveEffectByName("Islands") |
| 55 | +effect = segmentEditorWidget.activeEffect() |
| 56 | +segmentEditorNode.SetSelectedSegmentID(ID) |
| 57 | +effect.setParameterDefault("Operation", "KEEP_LARGEST_ISLAND") |
| 58 | +effect.self().onApply()#apply separate |
| 59 | +segmentEditorWidget.setActiveEffectByName("Islands") |
| 60 | +effect = segmentEditorWidget.activeEffect() |
| 61 | +segmentEditorNode.SetSelectedSegmentID(ID) |
| 62 | +effect.setParameterDefault("Operation", "KEEP_LARGEST_ISLAND") |
| 63 | +effect.self().onApply()#apply separate |
| 64 | + |
| 65 | +#STEP 3 |
| 66 | +##########################################export segmentation to model####################################### |
| 67 | +segmentationNode = getNode(ID) |
| 68 | +shNode = slicer.mrmlScene.GetSubjectHierarchyNode() |
| 69 | +exportFolderItemId = shNode.CreateFolderItem(shNode.GetSceneItemID(), "Segments") |
| 70 | +slicer.modules.segmentations.logic().ExportAllSegmentsToModels(segmentationNode, exportFolderItemId) |
| 71 | + |
| 72 | +#STEP 4 |
| 73 | +##########################Generate a label-map volume from the segmentation of the skull########################## |
| 74 | +labelmapVolumeNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLLabelMapVolumeNode') |
| 75 | +slicer.modules.segmentations.logic().ExportVisibleSegmentsToLabelmapNode(segmentationNode, labelmapVolumeNode, masterVolumeNode) |
| 76 | +slicer.util.saveNode(labelmapVolumeNode, "c:/tmp/BodyComposition-label.nrrd") |
| 77 | + |
| 78 | +#STEP 5 |
| 79 | +###########################Create a medial surface from the label-map using the BinaryThinningImageFilter################## |
| 80 | +##note this is computationally expensive step- many take many minutes to complete ##### |
| 81 | +import SampleData |
| 82 | +import SimpleITK as sitk |
| 83 | +import sitkUtils |
| 84 | +# Get input volume node |
| 85 | +inputVolumeNode_label = slicer.util.getNode('LabelMapVolume')#the string here is what you have the labelmap named as |
| 86 | +# Create new volume node for output |
| 87 | +binarythinning_outputVolumeNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode", "binarythinning_filter")#second string is name you give it |
| 88 | +# Run filter |
| 89 | +inputImage = sitkUtils.PullVolumeFromSlicer(inputVolumeNode_label) |
| 90 | +filter = sitk.BinaryThinningImageFilter() |
| 91 | +outputImage = filter.Execute(inputImage) |
| 92 | +sitkUtils.PushVolumeToSlicer(outputImage, binarythinning_outputVolumeNode) |
| 93 | +# Show processing result |
| 94 | +#slicer.util.setSliceViewerLayers(background=binarythinning_outputVolumeNode) |
| 95 | + |
| 96 | +#STEP 6 |
| 97 | +#############################Create a distance map from the medial surface using the DanielssonDistanceMapImageFilter ###############3 |
| 98 | +### with options Input is Binary = Yes and Use Image Spacing=Yes |
| 99 | +# Get input volume node |
| 100 | +inputVolumeNode_label = slicer.util.getNode("binarythinning_filter")#the string here "binarythinning_filter" is what you are naming the labelmap |
| 101 | +# Create new volume node for output |
| 102 | +DanielssonDistanceMap_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode", "DanielssonDistance")#second string is name you give it |
| 103 | + |
| 104 | +# Run filter |
| 105 | +inputImage = sitkUtils.PullVolumeFromSlicer(inputVolumeNode_label) |
| 106 | +filter = sitk.DanielssonDistanceMapImageFilter() |
| 107 | +filter.InputIsBinaryOn()#input is binary, |
| 108 | +filter.SetUseImageSpacing(True)#use image spaceing |
| 109 | +outputImage = filter.Execute(inputImage) |
| 110 | +sitkUtils.PushVolumeToSlicer(outputImage, DanielssonDistanceMap_node) |
| 111 | + |
| 112 | +#double the medial thickness to get a volume with total thickness, not medial thickness |
| 113 | +volumeNode = getNode('DanielssonDistance') |
| 114 | +a = arrayFromVolume(volumeNode) |
| 115 | +# Increase scalar values by two |
| 116 | +a[:] = a * 2.0 |
| 117 | +arrayFromVolumeModified(volumeNode) |
| 118 | +#type a again and see that the array values have doubled |
| 119 | + |
| 120 | + |
| 121 | +#STEP 7 |
| 122 | +#################### create the model with the colored thickness########################## |
| 123 | +#this function maps the thickness values onto the model. It essentially connects the |
| 124 | +#relevant parameters (i.e., input volume, model) into the existing C++ function |
| 125 | +#for mor information see https://slicer.readthedocs.io/en/latest/developer_guide/python_faq.html?highlight=CLI#how-to-run-a-cli-module-from-python |
| 126 | + |
| 127 | +def createprobevoltomodel(inputVolumeNode): |
| 128 | + """create texture on model using CLI module""" |
| 129 | + # Set parameters |
| 130 | + parameters = {} |
| 131 | + parameters["InputVolume"] = inputVolumeNode#input volume [image] |
| 132 | + parameters["InputModel"] = slicer.util.getFirstNodeByName(ID, className="vtkMRMLModelNode") |
| 133 | + parameters["OutputArrayName"] = "thickness_model" |
| 134 | + outputModelNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLModelNode","thickness") |
| 135 | + parameters["OutputModel"] = outputModelNode |
| 136 | + |
| 137 | + #parameters of CLU module |
| 138 | + #Group: IO |
| 139 | + # InputVolume [image]: Input volume |
| 140 | + # InputModel [geometry]: Input model |
| 141 | + # OutputModel [geometry]: Output model |
| 142 | + # OutputArrayName [string]: Output array name |
| 143 | + # Execute |
| 144 | + grayMaker = slicer.modules.probevolumewithmodel |
| 145 | + cliNode = slicer.cli.runSync(grayMaker, None, parameters) |
| 146 | + # Process results |
| 147 | + if cliNode.GetStatus() & cliNode.ErrorsMask: |
| 148 | + # error |
| 149 | + errorText = cliNode.GetErrorText() |
| 150 | + slicer.mrmlScene.RemoveNode(cliNode) |
| 151 | + raise ValueError("CLI execution failed: " + errorText) |
| 152 | + # success |
| 153 | + slicer.mrmlScene.RemoveNode(cliNode) |
| 154 | + return outputModelNode |
| 155 | + |
| 156 | +volumeNode = getNode("DanielssonDistance")#get the appropriate volume node |
| 157 | +createprobevoltomodel(volumeNode)#run the function and generate the model |
| 158 | + |
| 159 | + |
| 160 | +##STEP 8 |
| 161 | +#####Create the colormap label |
| 162 | +thicknessMapNode = getNode("thickness")#call the model with the thickness painted on it |
| 163 | + |
| 164 | +# Create color node |
| 165 | +colorNode = slicer.mrmlScene.CreateNodeByClass("vtkMRMLProceduralColorNode") |
| 166 | +colorNode.UnRegister(None) # to prevent memory leaks |
| 167 | +colorNode.SetName(slicer.mrmlScene.GenerateUniqueName("MedialThicknessMap")) |
| 168 | +#colorNode.SetAttribute("Category", "MedialThicknessModule") |
| 169 | +# The color node is a procedural color node, which is saved using a storage node. |
| 170 | +# Hidden nodes are not saved if they use a storage node, therefore the color node must be set to visible. |
| 171 | +colorNode.SetHideFromEditors(False) |
| 172 | +slicer.mrmlScene.AddNode(colorNode) |
| 173 | + |
| 174 | +#make |
| 175 | +import numpy as np |
| 176 | +scalars = vtk.util.numpy_support.vtk_to_numpy(thicknessMapNode.GetPolyData().GetPointData().GetScalars()) |
| 177 | +maxThickness = np.max(scalars) |
| 178 | +#set threshrange |
| 179 | +thresh = mythresh/maxThickness# this will be used to compute wheree the thickness threshold will occur on the label (as a proportion of max thickness) |
| 180 | + |
| 181 | +# Colormap - threshold |
| 182 | +colorMap = colorNode.GetColorTransferFunction() |
| 183 | +colorMap.RemoveAllPoints()#specify points for color map |
| 184 | +colorMap.AddRGBPoint(0.0, 0.0, 0.0, 1.0) |
| 185 | +colorMap.AddRGBPoint(thresh, 0.0, 0.0, 1.0)#put one color for all values below threshold |
| 186 | +colorMap.AddRGBPoint(thresh+0.01, 1.0, 0.0, 0.0)#put one color for all values below threshold |
| 187 | +colorMap.AddRGBPoint(1, 1.0, 0.0, 0.0)# |
| 188 | + |
| 189 | +# Display color legend |
| 190 | +thicknessMapNode.GetDisplayNode().SetAndObserveColorNodeID(colorNode.GetID()) |
| 191 | +##NOTE - if you want other default color schemes, just replace this line with the preceding. Here is the 'rainbow' theme: |
| 192 | +#thicknessMapNode.GetDisplayNode().SetAndObserveColorNodeID("vtkMRMLColorTableNodeRainbow") |
| 193 | + |
| 194 | +colorLegendDisplayNode = slicer.modules.colors.logic().AddDefaultColorLegendDisplayNode(thicknessMapNode) |
| 195 | +colorLegendDisplayNode.SetTitleText("Thickness (mm)") |
| 196 | +colorLegendDisplayNode.SetLabelFormat("%4.1f mm") |
| 197 | + |
| 198 | + |
0 commit comments