Volume Cropping in Python Script

Dear all,
I just checked Python_scripting and ScriptRepository pages and couldn’t find an example for cropping an image, I wrote this simple example. It runs without problem inside a module or extension but when I run it in an external python script it also runs and produces the correct output cropped image but there are error messages I don’t understand. Here is the mentioned complete script code:

 from __main__ import vtk, qt, ctk, slicer
 import sitkUtils

 print("                   myScript ")

 inputVolumeFnm  = ./myVolume.nrrd"  
 croppedImageFnm = "./myCroppedVolume.nrrd"

 [success, inputVolume] = slicer.util.loadVolume(inputVolumeFnm, returnNode=True)          
 inputImage = sitkUtils.PullVolumeFromSlicer(inputVolume.GetID())

 cropper = sitkUtils.sitk.CropImageFilter()
 croppingBounds = [[178, 210, 67],[227, 195, 34]]  
 croppedImage = cropper.Execute(inputImage, croppingBounds[0], croppingBounds[1]) 
 croppedNode = sitkUtils.PushVolumeToSlicer(croppedImage, None,  inputVolume.GetName() , 
 'vtkMRMLScalarVolumeNode' )

 properties = {}
 properties["fileType"] = ".nrrd"
 slicer.util.saveNode(croppedNode, croppedImageFnm, properties)

 print(" All tasks are done!  ") 

 exit()

Despite the correct output cropped image, I get this output in the terminal:

~/sw/Slicer-4.8.1/Slicer --no-main-window --python-script ./myScript.py
Qt: Cannot set locale modifiers: 
Number of registered modules: 150 
Number of instantiated modules: 150 
Number of loaded modules: 150 
               myScript 
Loaded volume from file: ./myVolume.nrrd. Dimensions: 485x485x121. Number of components: 1. Pixel type: short.
"Volume" Reader has successfully read the file     "./myVolume.nrrd" "[0.56s]" 
vtkITKArchetypeImageSeriesReader::ExecuteInformation: Archetype file ./_DUMMY_DOES_NOT_EXIST__ does not exist.
Algorithm vtkITKArchetypeDiffusionTensorImageReaderFile(0x4a977b0) returned failure for request: 
vtkInformation (0x4a9aea0)
 Debug: Off
 Modified Time: 827149
 Reference Count: 1
 Registered Events: (none)
 Request: REQUEST_INFORMATION
 ALGORITHM_AFTER_FORWARD: 1
 FORWARD_DIRECTION: 0
vtkITKArchetypeImageSeriesReader::ExecuteInformation: Archetype file ./_DUMMY_DOES_NOT_EXIST__ does not exist.
 Algorithm vtkITKArchetypeImageSeriesVectorReaderSeries(0x48edfb0) returned failure for request: vtkInformation (0x4a978f0)
  Debug: Off
  Modified Time: 828705
  Reference Count: 1
  Registered Events: (none)
 Request: REQUEST_INFORMATION
  ALGORITHM_AFTER_FORWARD: 1
  FORWARD_DIRECTION: 0
 vtkITKArchetypeImageSeriesReader::ExecuteInformation: Archetype file /_DUMMY_DOES_NOT_EXIST__ 
 does not exist.
 Algorithm vtkITKArchetypeImageSeriesScalarReader(0x4a6be60) returned failure for request: vtkInformation (0x4a8c600)
  Debug: Off
  Modified Time: 829510
  Reference Count: 1
  Registered Events: (none)
  Request: REQUEST_INFORMATION
  ALGORITHM_AFTER_FORWARD: 1
  FORWARD_DIRECTION: 0
  ReadData: This is not a nrrd file
  ReadData: Cannot read file as a volume of type DiffusionTensorVolume[fullName 
         /_DUMMY_DOES_NOT_EXIST__]
Number of files listed in the node = 0.
File reader says it was able to read 0 files.
File reader used the archetype file name of //_DUMMY_DOES_NOT_EXIST__ []
  FileNotFoundError
  ReadData: This is not a nrrd file
  ReadData: Failed to instantiate a file reader
  ReadData: Cannot read file as a volume of type Volume[fullName = /_DUMMY_DOES_NOT_EXIST__]
Number of files listed in the node = 0.
File reader says it was able to read 0 files.
File reader used the archetype file name of /_DUMMY_DOES_NOT_EXIST__ []
   FileNotFoundError
  Generic Warning: In /home/kitware/Dashboards/Package/Slicer- 
 481/Libs/MRML/Core/vtkDataFileFormatHelper.cxx, line 237
 vtkDataFileFormatHelper::GetFileExtensionFromFormatString: please update deprecated extension-only 
 format specifier to 'File format name (.ext)' format! Current format string: .nrrd
 All tasks are done!
1 Like

This warning is generated by a workaround in sitkUtils. You can safely ignore the warnings.

@blowekamp Do you think we still need this workaround? On Windows, it does not seem to be necessary - for example, this code works:

import SimpleITK as sitk

# Pull volume from Slicer
inputVolumeNode = getNode('MRHead')
inputVolumeNodeSceneAddress = inputVolumeNode.GetScene().GetAddressAsString("").replace('Addr=','')
inputVolumeNodeFullITKAddress = 'slicer:{0}#{1}'.format(inputVolumeNodeSceneAddress,inputVolumeNode.GetID())
sitkimage = sitk.ReadImage(inputVolumeNodeFullITKAddress)

# Process
filter = sitk.SignedMaurerDistanceMapImageFilter()
outputImage = filter.Execute(sitkimage)

# Push volume to Slicer
outputVolumeNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode")
outputVolumeNode.CreateDefaultDisplayNodes()
outputVolumeNodeSceneAddress = outputVolumeNode.GetScene().GetAddressAsString("").replace('Addr=','')
outputVolumeNodeFullITKAddress = 'slicer:{0}#{1}'.format(inputVolumeNodeSceneAddress,outputVolumeNode.GetID())
sitk.WriteImage(outputImage, outputVolumeNodeFullITKAddress)
slicer.util.setSliceViewerLayers(background = outputVolumeNode)
2 Likes

That code appears to be added by yourself:

The ITK factory registration is a complicated process and likely has changed a bit since I last looked at it in Slicer so I don’t know.

With SimpleITK 1.1 we have added a method to get a list of the registered ITK ImageIO. This could detect if this need to happen or not.
https://itk.org/SimpleITKDoxygen110/html/classitk_1_1simple_1_1ImageReaderBase.html#a6d5066467465ebff4b2bb895e7ed59d0

I’ve tried to access std::vector<std::string> itk::simple::ImageReaderBase::GetRegisteredImageIOs() from Python but it did not seem to be available. Do you know if std::vector<std::string> types are supposed to be wrapped for Python in SimpleITK? (I know that VTK wrappers cannot deal with such types)

I just reorganized the code, the code for ensuring registration was added by @hjmjohnson.

Currently, Slicer is using 1.0.1-ish there is a pull request to update to 1.1.0..

If you get SimpleITK with pip it will be there.

Yes the return type is wrapped and tested in SimpleITK.


In [1]: import SimpleITK as sitk

In [2]: print sitk.Version()
SimpleITK Version: 1.1.0 (ITK 4.13)
Compiled: Mar 22 2018 21:14:44


In [3]: print sitk.ImageFileReader().GetRegisteredImageIOs()
('BMPImageIO', 'BioRadImageIO', 'Bruker2dseqImageIO', 'GDCMImageIO', 'GE4ImageIO', 'GE5ImageIO', 'GiplImageIO', 'HDF5ImageIO', 'JPEGImageIO', 'LSMImageIO', 'MINCImageIO', 'MRCImageIO', 'MetaImageIO', 'NiftiImageIO', 'NrrdImageIO', 'PNGImageIO', 'StimulateImageIO', 'TIFFImageIO', 'VTKImageIO')

Great! I’ll test Slicer with SimpleITK 1.1.0 tomorrow and if everything works then I’ll update the nightly build to use that version.

1 Like

Hello @brhoom, I saw your post here and thought you might be familiar with cropping volume in Python. I am wondering if you know of a way to crop a volume using ROI in 3D slicer in Python, modify the cropped volume and put it back into the original volume (at the same position). Essentially, what I am trying to do is like this

Drag some ROI in 3D slicer
croppedVolumeNode = some_internal_cropping_function(volumeNode, roiNode)
modify the cropped volume
volumeNode.update(croppedVolumeNode) # to put it back to the original location

Thanks!

Sorry for the late reply. I do the cropping using only one point, the center of the cropping. So to get the croppingBounds in my code above from a defined point:

    spacing = inputVolume.GetSpacing()
    imgData = inputVolume.GetImageData()
    dimensions = imgData.GetDimensions() 
    size = [lengthX,lengthY,lengthZ] # defined volume dimensions    
    # Calculate lower and upper cropping bounds
    lower = [0,0,0]
    for i in range(0,3):
        lower[i] = point[i] - int((size[i]/spacing[i])/2)
    upper = [0,0,0]
    for i in range(0,3):
        upper[i] = dimensions[i] - (point[i]+int((size[i]/spacing[i])/2))
    croppingBounds = [lower,upper]

I think in your case, I will do the following :slight_smile:

  • create a copy of the original volume.
  • do the cropping
  • process the cropped volume.
  • convert the copied volume and the cropped volume to np arrays and replace the voxels in the first volume with the cropped one.

Probably there is a better way. I am not the Slicer expert here. Also, I think it is not a good idea to ask different question in different post. I would open a new post with the new question to keep the forum organized.

Hope this helps.

Yeah I posted a separate question in the forum. Thank you for your reply anyway!