Slicer CPPYY example

guide.py

# Copyright Steve Pieper
# https://github.com/Slicer/Slicer/blob/master/License.txt
# Demo video: https://youtu.be/xcQKj4yp2nw

"""

exec(open('/Users/pieper/slicer/latest/SlicerTMSGuide/guide.py', 'r').read())

"""

import cppyy
import numpy
import os
import time

def arrayFromModelPointDataModified(modelNode, arrayName):
  """Indicate that modification of a numpy array returned by :py:meth:`arrayFromModelPoints` has been completed."""
  arrayVtk = modelNode.GetPolyData().GetPointData().GetArray(arrayName)
  arrayVtk.Modified()

slicer.util.arrayFromModelPointDataModified =  arrayFromModelPointDataModified

slicer.mrmlScene.Clear(0)

path = '/Users/pieper/slicer/latest/SlicerTMSGuide'
pialPath = os.path.join(path, 'data/lh.pial')
thicknessPath = os.path.join(path, 'data/lh.thickness')

pial = slicer.util.loadModel(pialPath)
slicer.util.loadScalarOverlay(thicknessPath, pial.GetID())

pial.GetDisplayNode().SetScalarRangeFlag(0)
pial.GetDisplayNode().SetScalarRange(0,50)
pial.GetDisplayNode().SetAndObserveColorNodeID('vtkMRMLColorTableNodeFileCividis.txt')

ras = [-65, 10, 26]
markupsLogic = slicer.vtkSlicerMarkupsLogic()
markupsLogic.SetMRMLScene(slicer.mrmlScene)
markupsLogic.AddFiducial(*ras)
markupsLogic.JumpSlicesToLocation(*ras, False)
markupFiducials = slicer.util.getNode('F')
markupFiducials.GetDisplayNode().SetTextScale(0)
markupFiducials.GetDisplayNode().SetGlyphScale(10)
markupFiducials.GetDisplayNode().SetSelectedColor([1,1,0])
markupFiducials.GetDisplayNode().SetSnapMode(slicer.vtkMRMLMarkupsDisplayNode.SnapModeUnconstrained)

def updateOverlay(caller, event):
    pointsArray = slicer.util.arrayFromModelPoints(pial)
    scalarArray = slicer.util.arrayFromModelPointData(pial, 'lh.thickness')
    pointArray = slicer.util.arrayFromMarkupsControlPoints(markupFiducials)
    func1d = lambda p: numpy.linalg.norm(pointArray - p)
    scalarArray[:] = numpy.apply_along_axis(func1d, 1, pointsArray)
    slicer.util.arrayFromModelPointDataModified(pial, 'lh.thickness')

def compileAndImportCPP(cppSource):
    namespaceTag = "cppyy_"+str(time.time()).replace(".", "_")
    cppSource = "namespace " + namespaceTag + " {" + cppSource + " }"
    cppyy.cppdef(cppSource)
    exec("from cppyy.gbl import " + namespaceTag)
    namespace = eval("cppyy.gbl." + namespaceTag)
    return namespace

cppyy.include("math.h")
fillWithDistancesSource = """
    void fillWithDistances(void *pointsArrayData, void *pointData, void *scalarArrayData, unsigned int size) {
        float *pointsArray = (float *)pointsArrayData;
        double *point = (double *)pointData;
        float *scalarArray = (float *)scalarArrayData;
        unsigned int offset;
        for (offset = 0; offset < size; offset++) {
            scalarArray[offset] = sqrt ( 
              (pointsArray[offset*3 + 0] - point[0]) * (pointsArray[offset*3 + 0] - point[0]) +
              (pointsArray[offset*3 + 1] - point[1]) * (pointsArray[offset*3 + 1] - point[1]) +
              (pointsArray[offset*3 + 2] - point[2]) * (pointsArray[offset*3 + 2] - point[2]) );
        }
    }
"""
fillWithDistances = compileAndImportCPP(fillWithDistancesSource).fillWithDistances

def updateOverlayCpp(caller, event):
    pointsArray = slicer.util.arrayFromModelPoints(pial)
    scalarArray = slicer.util.arrayFromModelPointData(pial, 'lh.thickness')
    pointArray = slicer.util.arrayFromMarkupsControlPoints(markupFiducials)
    fillWithDistances(pointsArray.data, pointArray.data, scalarArray.data, scalarArray.size)
    slicer.util.arrayFromModelPointDataModified(pial, 'lh.thickness')

pointsArray = slicer.util.arrayFromModelPoints(pial)
scalarArray = slicer.util.arrayFromModelPointData(pial, 'lh.thickness')
pointArray = slicer.util.arrayFromMarkupsControlPoints(markupFiducials)
updateOverlayCpp(markupFiducials, slicer.vtkMRMLMarkupsNode.PointModifiedEvent)

ras = [50, 0, -70]
labelFiducial = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode")
markupsLogic.SetActiveListID(labelFiducial)
markupsLogic.AddFiducial(*ras)
markupsLogic.JumpSlicesToLocation(*ras, False)
labelFiducial.GetDisplayNode().SetTextScale(10)
labelFiducial.GetDisplayNode().SetGlyphScale(0)
labelFiducial.GetDisplayNode().SetSelectedColor([1,1,0])
labelFiducial.GetDisplayNode().SetSnapMode(slicer.vtkMRMLMarkupsDisplayNode.SnapModeUnconstrained)

#mode = "python"
mode = "cppyy"
if mode == "python":
    labelFiducial.SetNthControlPointLabel(0, "Python")
    markupFiducials.AddObserver(slicer.vtkMRMLMarkupsNode.PointModifiedEvent, updateOverlay)
else:
    labelFiducial.SetNthControlPointLabel(0, "Python + CPPYY")
    markupFiducials.AddObserver(slicer.vtkMRMLMarkupsNode.PointModifiedEvent, updateOverlayCpp)

Very nice! I would recommend to add pip_install("cppyy") to the script and download all input data automatically to ensure that it works just by copy-pasting the code into the Python interactor.

Very nice!

Thanks!

I would recommend to add pip_install("cppyy") to the script and download all input data automatically to ensure that it works just by copy-pasting the code into the Python interactor.

Yes, good points - a caveat here is that cppyy is still experimental and only works on a developer build for now since the install step requires compilation. So this is not really ready to promote to the whole community, but it's close.

download all input data automatically

We are working an updated neuroimaging dataset but I don't have it released in sampledata yet. Of course this demo could work with any model.

Maybe the best thing is just to strip down the scalar array part to a small example and put it in the script repository for reference.

Installing on a mac build currently requires setting this variable so that the pip install succeeds in building the prereqs: SDKROOT=/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk ./Slicer-build/Slicer