Python interactor

The Python interactor is a very useful utility to access run-time any 3DSlicer and SlicerAstro class: MRML nodes, Widgets and Logics (N.B. the data of the fits file are stored in the MRML AstroVolume class). Python scripts can also be run by typing in a terminal window:

Slicer --python-script example.py

or by copying and pasting the code in the python console of 3DSlicer:

See also:

• https://www.slicer.org/slicerWiki/index.php/Documentation/Nightly/Developers/Python_scripting

Note:

• Slicer (and therefore also SlicerAstro) can be used as a fully interactive jupyter notebook, see link.

Install Additional Python packages

3DSlicer and SlicerAstro binaries deliver their own Python environment. The Python version is 3.6.7 and it incorporates the additional packages numpy, astropy, scipy. It is possible to install additional 'pure' python wheels. Instructions are at the following link.

3DSlicer Python scripts repository

• https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository

SlicerAstro Python scripts repository

as an example a number of Python scripts are provided below for a variety of applications:

• SmoothRenderSaveVideos.py: This script requires ffmpeg which for license incompatibility can not be included in 3DSlicer binaries. It is therefore necessary to install the ffmpeg binary manually when it is not present in the system, and select it in the advanced widget of the screen capture (this last step is required only once).

  slicer.util.loadVolume("/full_path/galaxy.fits",{"center":True})
  mw = slicer.util.mainWindow()
  ms = mw.moduleSelector()
  # Smooth the datacube
  ms.selectModule('AstroSmoothing')
  smowidget = slicer.modules.astrosmoothing.widgetRepresentation()
  smomrmlpara = smowidget.mrmlAstroSmoothingParametersNode()
  smomrmlpara.SetMode("Manual")
  # only if a GPU is present:
  #smomrmlpara.SetHardware(1)
  smowidget.onApply()
  ms.selectModule('AstroVolume')
  # Setting maximum quality for the rendering
  astrovolumewidget = slicer.modules.astrovolume.widgetRepresentation()
  astrovolumewidget.onCurrentQualityControlChanged(1)
  volumes = slicer.mrmlScene.GetNodesByClass("vtkMRMLAstroVolumeNode")
  volumefiltered = volumes.GetItemAsObject(1)
  smomrmlpara.SetInputVolumeNodeID(volumefiltered.GetID())
  astrovolumewidget.onCurrentQualityControlChanged(1)
  # Create videos 
  ms.selectModule('ScreenCapture')
  screencapturewidget = slicer.modules.screencapture.widgetRepresentation()
  instance = screencapturewidget.self()
  viewNode = slicer.util.getNode('vtkMRMLViewNode1')
  instance.viewNodeSelector.setCurrentNode(viewNode)
  instance.numberOfStepsSliderWidget.setValue(360)
  instance.videoExportCheckBox.setChecked(1)
  instance.videoFormatWidget.setCurrentIndex(1)
  instance.videoFileNameWidget.setText("WEIN069.mp4")
  instance.videoLengthSliderWidget.setValue(6)
  instance.onCaptureButton()
  viewNode = slicer.util.getNode('vtkMRMLViewNode2')
  instance.viewNodeSelector.setCurrentNode(viewNode)
  instance.viewNodeSelector.setCurrentNode(viewNode)
  instance.numberOfStepsSliderWidget.setValue(360)
  instance.videoExportCheckBox.setChecked(1)
  instance.videoFormatWidget.setCurrentIndex(1)
  instance.videoFileNameWidget.setText("WEIN069_smoothed.mp4")
  instance.videoLengthSliderWidget.setValue(6)
  instance.onCaptureButton()

• Create(alike)MomentMap.py:

  # Get the data-cube volume
  slicer.util.loadVolume("/full_path/galaxy.fits",{"center":True})
  volumeNode = getNode('galaxy')
  datacube = volumeNode.GetImageData()
  
  # Get dimensions
  N1 = int(volumeNode.GetAttribute("SlicerAstro.NAXIS1"))
  N2 = int(volumeNode.GetAttribute("SlicerAstro.NAXIS2"))
  RMS = float(volumeNode.GetAttribute("SlicerAstro.DisplayThreshold"))
  
  # Create an empty 2-D image
  imageSize = [N1, N2, 1]
  imageSpacing = [1.0, 1.0, 1.0]
  voxelType = vtk.VTK_FLOAT
  imageDataTemp = vtk.vtkImageData()
  imageDataTemp.SetDimensions(imageSize)
  imageDataTemp.SetSpacing(imageSpacing)
  imageDataTemp.AllocateScalars(voxelType, 1)
  
  extentTemp = imageDataTemp.GetExtent()
  for i in range(extentTemp[0], extentTemp[1]+1):
     for j in range(extentTemp[2], extentTemp[3]+1):
        for k in range(extentTemp[4], extentTemp[5]+1):
           imageDataTemp.SetScalarComponentFromFloat(i,j,k,0,0.)
  
  # calculate moment map
  imageData = volumeNode.GetImageData()
  extent = imageData.GetExtent()
  for i in range(extent[0], extent[1]+1):
     for j in range(extent[2], extent[3]+1):
        sum = 0.
        for k in range(extent[4], extent[5]+1):
           value = imageData.GetScalarComponentAsFloat(i,j,k,0)
           if value > 3 * RMS:
              sum += value
        imageDataTemp.SetScalarComponentFromFloat(i,j,0,0,sum)
  
  imageDataTemp.Modified()
  point = imageDataTemp.GetPointData()
  array = point.GetArray("ImageScalars")
  point.Modified()
  array.Modified()
  array.GetValueRange()
  
  # create Astro Volume for the moment map
  astroVolumeLogic = slicer.modules.astrovolume.logic()
  volumeNodeMomentMap = astroVolumeLogic.CloneVolume(slicer.mrmlScene, volumeNode, 'MomentMap')
  
  # modify fits attributes
  volumeNodeMomentMap.SetAttribute("SlicerAstro.NAXIS", "2")
  volumeNodeMomentMap.RemoveAttribute("SlicerAstro.NAXIS3")
  volumeNodeMomentMap.RemoveAttribute("SlicerAstro.CROTA3")
  volumeNodeMomentMap.RemoveAttribute("SlicerAstro.CRPIX3")
  volumeNodeMomentMap.RemoveAttribute("SlicerAstro.CRVAL3")
  volumeNodeMomentMap.RemoveAttribute("SlicerAstro.CTYPE3")
  volumeNodeMomentMap.RemoveAttribute("SlicerAstro.CUNIT3")
  volumeNodeMomentMap.RemoveAttribute("SlicerAstro.DTYPE3")
  volumeNodeMomentMap.RemoveAttribute("SlicerAstro.DRVAL3")
  volumeNodeMomentMap.RemoveAttribute("SlicerAstro.DUNIT3")
  
  # copy 2-D image into the Astro Volume object
  volumeNodeMomentMap.SetAndObserveImageData(imageDataTemp)
  volumeNodeMomentMap.UpdateDisplayThresholdAttributes()
  volumeNodeMomentMap.UpdateRangeAttributes()
  
  # change colorMap of the 2-D image
  displayNode = volumeNodeMomentMap.GetDisplayNode()
  displayNode.SetAndObserveColorNodeID('vtkMRMLColorTableNodeRainbow')
  

• AccessDataAsNumpy.py:

  def astroArrayFromVolume(volumeNode):
     """Return a voxel array from the volume node as a numpy array.
     Voxel values are not copied. Voxel values in the volume node can be modified
     by changing values in the numpy array.
     After all modifications have been completed, call volumeNode.Modified().
  
     .. warning:: the memory area of the returned array is managed by VTK, therefore
     values in the array may be changed, but the array must not be reallocated
     (change of array size, shallow-copy content from other array most likely will cause
     the application to crash). To allow arbitrary numpy operations on a volume array:
  
     1. Make a deep-copy of the returned VTK-managed array using :func:`numpy.copy`.
     2. Perform any computations using the copied array.
     3. Write results back to the image data using :py:meth:`updateVolumeFromArray`.
     """
  scalarTypes = ['vtkMRMLAstroVolumeNode', 'vtkMRMLAstroLabelMapVolumeNode']
  vimage = volumeNode.GetImageData()
  nshape = tuple(reversed(volumeNode.GetImageData().GetDimensions()))
  narray = None
  if volumeNode.GetClassName() in scalarTypes:
     narray = vtk.util.numpy_support.vtk_to_numpy(vimage.GetPointData().GetScalars()).reshape(nshape)       
  return narray
  
  

  def astroUpdateVolumeFromArray(volumeNode, narray):
     """Sets voxels of a volume node from a numpy array.
     Voxels values are deep-copied, therefore if the numpy array
     is modified after calling this method, voxel values in the volume node will not change.
     Dimensions and data size of the source numpy array do not have to match the current
     content of the volume node.
     """
  
  vshape = tuple(reversed(narray.shape))
  if len(vshape) == 1:
     narray = numpy.expand_dims(narray, axis=0)
     narray = numpy.expand_dims(narray, axis=0)
     vshape = tuple(reversed(narray.shape))
     vcomponents = 1
  elif len(vshape) == 2:
     narray = numpy.expand_dims(narray, axis=0)
     vshape = tuple(reversed(narray.shape))
     vcomponents = 1
  elif len(vshape) == 3:
     vcomponents = 1
  else:
     raise RuntimeError("Unsupported numpy array shape: "+str(narray.shape))
  
  vimage = volumeNode.GetImageData()
  if vimage is None:
     vimage = vtk.vtkImageData()
     volumeNode.SetAndObserveImageData(vimage)
  
  vtype = vtk.util.numpy_support.get_vtk_array_type(narray.dtype)
  vimage.SetDimensions(vshape)
  vimage.AllocateScalars(vtype, vcomponents)
  
  narrayTarget = astroArrayFromVolume(volumeNode)
  narrayTarget[:] = narray
  
  # Notify the application that image data is changed
  # (same notifications as in vtkMRMLVolumeNode.SetImageDataConnection)
  imageData = volumeNode.GetImageData()
  pointData = imageData.GetPointData() if imageData else None
  if pointData is not None:
     if pointData.GetScalars() is not None:
        pointData.GetScalars().Modified()
  
  volumeNode.StorableModified()
  volumeNode.Modified()
  volumeNode.InvokeEvent(slicer.vtkMRMLVolumeNode.ImageDataModifiedEvent, volumeNode)
  

• maskVolumeWithLabelVolume.py:

  # Get the data-cube volumes
  volumeNode = getNode('Name')
  volumeLabelNode = getNode('LabelName')
  
  # masking
  imageData = volumeNode.GetImageData()
  imageDataLabel = volumeLabelNode.GetImageData()
  extent = imageData.GetExtent()
  for i in range(extent[0], extent[1]+1):
     for j in range(extent[2], extent[3]+1):
        for k in range(extent[4], extent[5]+1):
           if imageDataLabel.GetScalarComponentAsFloat(i,j,k,0) < 1 :
              imageData.SetScalarComponentFromFloat(i,j,k,0, 0.)
  
  # update volume keywords
  volumeNode.UpdateDisplayThresholdAttributes()
  volumeNode.UpdateRangeAttributes()
  

• accessSlicerAstrovtkFits.py:

  import vtkFitsPython as vtkFits
  vtkFits.vtkFITSReader()
  dir(vtkFits)
  

• loadAndVisualizeSofiaCatalog.py:

  import numpy as np
  
  # load the data
  slicer.util.loadVolume("fullpath/sofiatestcube.fits",{"center":True})
  volumeNode = getNode("sofiatestcube")
  volumeDisplayNode = volumeNode.GetDisplayNode()
  IJKToRASMatrix = vtk.vtkMatrix4x4()
  volumeNode.GetIJKToRASMatrix(IJKToRASMatrix)
  
  # create markups node
  markupsNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode")
  markupsNode.CreateDefaultDisplayNodes()
  markupsDisplayNode = markupsNode.GetDisplayNode()
  markupsDisplayNode.SetGlyphScale(3)
  markupsDisplayNode.SetTextScale(3)
  markupsDisplayNode.SliceProjectionOn()
  markupsDisplayNode.SliceProjectionUseFiducialColorOff()
  markupsDisplayNode.SetColor(0.4, 1., 1.)
  markupsDisplayNode.SetSelectedColor(0.4, 1., 1.)
  markupsDisplayNode.SetSliceProjectionColor(1., 0.733, 0.078)
  
  # load the Sofia catalog
  table = np.genfromtxt("fullpath/test4slicerastro_cat.ascii", dtype='str')
  shape = table.shape
   
  # In catalogue files produced by Sofia version 1.2 the column indexes for  
  # ra, dec and vopt are respectively: 33, 34, 35.
  # It is assumed that the epoch of the catalog is the same as the epoch of the  
  # data on which the catalog is over-plotted (volumeNode)    
  for ii in range(shape[0]):
     world = [(float)(table[ii][32]), (float) (table[ii][33]), (float) (table[ii][34])]
     ijk = [0., 0., 0.]
     volumeDisplayNode.GetIJKSpace(world, ijk)
     position_Ijk = [ijk[0], ijk[1], ijk[2], 1]
     ras = IJKToRASMatrix.MultiplyPoint(position_Ijk)
     markupsNode.AddFiducial(ras[0], ras[1], ras[2])
     markupsNode.SetNthFiducialLabel(ii, table[ii][36])
     markupsNode.SetNthMarkupLocked(ii, 1)