Note
Go to the end to download the full example code.
Cupy vs SitkTransforms#
author: OpenTPS team
This example shows how to use the OpenTPS transforms with cupy and compare the results with the SimpleITK transforms.
running time: ~ 5 minutes
Setting up the environment in google collab#
You will first need to change the execution parameters in your Colab environment (bottom left) from CPU to GPU so that you can run the example.
import sys
if "google.colab" in sys.modules:
from IPython import get_ipython
get_ipython().system('git clone https://gitlab.com/openmcsquare/opentps.git')
get_ipython().system('pip install ./opentps')
get_ipython().system('pip install cupy-cuda12x')
import opentps
imports
import copy
import matplotlib.pyplot as plt
import logging
import numpy as np
import os
import the needed opentps.core packages
from opentps.core.data.images import VectorField3D
from opentps.core.data.images import CTImage
from opentps.core.data.images import ROIMask
from opentps.core.data._transform3D import Transform3D
from opentps.core.processing.imageProcessing.imageTransform3D import rotateData, translateData
logger = logging.getLogger(__name__)
Output path#
output_path = os.path.join(os.getcwd(), 'Output', 'CupyVsSitk')
if not os.path.exists(output_path):
os.makedirs(output_path)
logger.info('Files will be stored in {}'.format(output_path))
Generate synthetic input immages#
imgSize = [40, 40, 40]
imgSpacing = [1, 1, 2]
objectBorder = [[21, 33], [int(imgSize[1]/4), 3*int(imgSize[1]/4)], [21, 34]]
translation = np.array([-10.22, 0, -14.56])
rotation = np.array([0, 30, 0])
rotCenter = 'imgCenter'
outputBox = 'same'
# interpOrder = 1
showImage = True
showField = True
showMask = True
fixed = CTImage()
fixed.spacing = np.array(imgSpacing)
fixed.imageArray = np.full(imgSize, -1000)
fixed.imageArray[objectBorder[0][0]: objectBorder[0][1],
objectBorder[1][0]: objectBorder[1][1],
objectBorder[2][0]: objectBorder[2][1]] = 100.0
y_slice = int(imgSize[1]/2)
pointList = [[objectBorder[0][0], y_slice, objectBorder[2][1]-1],
[objectBorder[0][1]-1, y_slice, objectBorder[2][0]],
[objectBorder[0][0]+1, y_slice, objectBorder[2][0]+1],
[objectBorder[0][0], y_slice, objectBorder[2][0]]]
fieldFixed = VectorField3D()
fieldFixed.imageArray = np.zeros((imgSize[0], imgSize[1], imgSize[2], 3))
fieldFixed.spacing = np.array(imgSpacing)
vectorList = [np.array([4, 6, 8]), np.array([0, 6, 8]), np.array([14, 6, 6]), np.array([6, 0, 0])]
for pointIdx in range(len(pointList)):
fieldFixed.imageArray[pointList[pointIdx][0], pointList[pointIdx][1], pointList[pointIdx][2]] = vectorList[
pointIdx]
maskFixed = ROIMask()
maskFixed.spacing = np.array(imgSpacing)
maskFixed.imageArray = np.zeros(imgSize).astype(bool)
maskFixed.imageArray[objectBorder[0][0]: objectBorder[0][1],
objectBorder[1][0]: objectBorder[1][1],
objectBorder[2][0]: objectBorder[2][1]] = True
this function is just to see the results
def showImagesAndFieldAndMask(fixed, movingCupy1, movingCupy3, movingSitk, fieldFixed, fieldMovingCupy1, fieldMovingCupy3, fieldMovingSitk,
maskFixed, maskMovingCupy1, maskMovingCupy3, maskMovingSitk, y_slice, figTitle, showImage=True,
showField=True, showMask=True, ):
compXFixed = fieldFixed.imageArray[:, y_slice, :, 0]
compZFixed = fieldFixed.imageArray[:, y_slice, :, 2]
compXMovingCupy1 = fieldMovingCupy1.imageArray[:, y_slice, :, 0]
compZMovingCupy1 = fieldMovingCupy1.imageArray[:, y_slice, :, 2]
compXMovingCupy3 = fieldMovingCupy3.imageArray[:, y_slice, :, 0]
compZMovingCupy3 = fieldMovingCupy3.imageArray[:, y_slice, :, 2]
compXMovingSITK = fieldMovingSitk.imageArray[:, y_slice, :, 0]
compZMovingSITK = fieldMovingSitk.imageArray[:, y_slice, :, 2]
fig, ax = plt.subplots(3, 4)
fig.suptitle(figTitle)
if showImage:
ax[0, 0].imshow(fixed.imageArray[:, y_slice, :])
ax[0, 1].imshow(movingCupy1.imageArray[:, y_slice, :])
ax[0, 2].imshow(movingCupy3.imageArray[:, y_slice, :])
ax[0, 3].imshow(movingSitk.imageArray[:, y_slice, :])
ax[1, 0].imshow(movingCupy1.imageArray[:, y_slice, :] - movingSitk.imageArray[:, y_slice, :])
ax[1, 0].set_xlabel('Img diff cupy1-sitk')
ax[1, 1].imshow(movingCupy1.imageArray[:, y_slice, :] - movingCupy3.imageArray[:, y_slice, :])
ax[1, 1].set_xlabel('Img diff cupy1-cupy3')
if showField:
ax[0, 0].quiver(compZFixed, compXFixed, alpha=0.5, color='red', angles='xy', scale_units='xy', scale=2, width=.010)
ax[0, 1].quiver(compZMovingCupy1, compXMovingCupy1, alpha=0.5, color='green', angles='xy', scale_units='xy', scale=2, width=.010)
ax[0, 2].quiver(compZMovingCupy3, compXMovingCupy3, alpha=0.5, color='green', angles='xy', scale_units='xy', scale=2, width=.010)
ax[0, 3].quiver(compZMovingSITK, compXMovingSITK, alpha=0.5, color='green', angles='xy', scale_units='xy', scale=2, width=.010)
ax[1, 2].quiver(compZMovingCupy1 - compZMovingSITK, compXMovingCupy1 - compXMovingSITK, alpha=0.5, color='green', angles='xy', scale_units='xy', scale=2, width=.010)
ax[1, 2].set_xlabel('Field diff cupy1-sitk')
ax[1, 3].quiver(compZMovingCupy1 - compZMovingCupy3, compXMovingCupy1 - compXMovingCupy3, alpha=0.5, color='green', angles='xy', scale_units='xy', scale=2, width=.010)
ax[1, 3].set_xlabel('Field diff cupy1-cupy3')
if showMask:
ax[0, 0].imshow(maskFixed.getBinaryContourMask(internalBorder=True).imageArray[:, y_slice, :], alpha=0.5, cmap='Reds')
ax[0, 1].imshow(maskMovingCupy1.getBinaryContourMask(internalBorder=True).imageArray[:, y_slice, :], alpha=0.5, cmap='Reds')
ax[0, 2].imshow(maskMovingCupy3.getBinaryContourMask(internalBorder=True).imageArray[:, y_slice, :], alpha=0.5, cmap='Reds')
ax[0, 3].imshow(maskMovingSitk.getBinaryContourMask(internalBorder=True).imageArray[:, y_slice, :], alpha=0.5, cmap='Reds')
ax[2, 0].imshow(maskMovingCupy1.getBinaryContourMask(internalBorder=True).imageArray[:, y_slice, :] ^ maskMovingSitk.getBinaryContourMask(internalBorder=True).imageArray[:, y_slice, :], alpha=0.5, cmap='Reds')
ax[2, 0].set_xlabel('Mask diff cupy1-sitk')
ax[2, 1].imshow(maskMovingCupy1.getBinaryContourMask(internalBorder=True).imageArray[:, y_slice,:] ^ maskMovingCupy3.getBinaryContourMask(internalBorder=True).imageArray[:, y_slice, :], alpha=0.5, cmap='Reds')
ax[2, 1].set_xlabel('Mask diff cupy1-cupy3')
ax[0, 0].set_title('Fixed')
ax[0, 0].set_xlabel(f"{fixed.origin}\n{fixed.spacing}\n{fixed.gridSize}")
ax[0, 1].set_title('Moving Cupy1')
ax[0, 1].set_xlabel(f"{movingCupy1.origin}\n{movingCupy1.spacing}\n{movingCupy1.gridSize}")
ax[0, 2].set_title('Moving Cupy3')
ax[0, 2].set_xlabel(f"{movingCupy3.origin}\n{movingCupy3.spacing}\n{movingCupy3.gridSize}")
ax[0, 3].set_title('Moving SITK')
ax[0, 3].set_xlabel(f"{movingSitk.origin}\n{movingSitk.spacing}\n{movingSitk.gridSize}")
plt.savefig(os.path.join(output_path, 'cupy_VS_sitk_Transforms.png'))
plt.show()
Test using a Transform3D#
print('-' * 40)
movingCupy1 = copy.deepcopy(fixed)
movingCupy3 = copy.deepcopy(fixed)
movingSitk = copy.deepcopy(fixed)
fieldMovingCupy1 = copy.deepcopy(fieldFixed)
fieldMovingCupy3 = copy.deepcopy(fieldFixed)
fieldMovingSitk = copy.deepcopy(fieldFixed)
maskMovingCupy1 = copy.deepcopy(maskFixed)
maskMovingCupy3 = copy.deepcopy(maskFixed)
maskMovingSitk = copy.deepcopy(maskFixed)
## Create a transform 3D
print('Create a transform 3D')
transform3D = Transform3D()
transform3D.initFromTranslationAndRotationVectors(transVec=translation, rotVec=rotation)
transform3D.setCenter(rotCenter)
print('Translation', transform3D.getTranslation())
print('Rotation', transform3D.getRotationAngles(inDegrees=True))
print('Moving with transform3D')
movingCupy1 = transform3D.deformData(movingCupy1, outputBox=outputBox, fillValue=-1000, tryGPU=True, interpOrder=1)
fieldMovingCupy1 = transform3D.deformData(fieldMovingCupy1, outputBox=outputBox, tryGPU=True, interpOrder=1)
maskMovingCupy1 = transform3D.deformData(maskMovingCupy1, outputBox=outputBox, tryGPU=True, interpOrder=1)
movingCupy3 = transform3D.deformData(movingCupy3, outputBox=outputBox, fillValue=-1000, tryGPU=True, interpOrder=3)
fieldMovingCupy3 = transform3D.deformData(fieldMovingCupy3, outputBox=outputBox, tryGPU=True, interpOrder=3)
maskMovingCupy3 = transform3D.deformData(maskMovingCupy3, outputBox=outputBox, tryGPU=True, interpOrder=3)
movingSitk = transform3D.deformData(movingSitk, outputBox=outputBox, fillValue=-1000)
fieldMovingSitk = transform3D.deformData(fieldMovingSitk, outputBox=outputBox)
maskMovingSitk = transform3D.deformData(maskMovingSitk, outputBox=outputBox)
showImagesAndFieldAndMask(fixed, movingCupy1, movingCupy3, movingSitk, fieldFixed, fieldMovingCupy1,
fieldMovingCupy3, fieldMovingSitk,
maskFixed, maskMovingCupy1, maskMovingCupy3, maskMovingSitk, y_slice,
figTitle='Test using a Transform3D', showField=showField, showImage=showImage,
showMask=showMask)
Test using translateData#
print('-' * 40)
movingCupy1 = copy.deepcopy(fixed)
movingCupy3 = copy.deepcopy(fixed)
movingSitk = copy.deepcopy(fixed)
fieldMovingCupy1 = copy.deepcopy(fieldFixed)
fieldMovingCupy3 = copy.deepcopy(fieldFixed)
fieldMovingSitk = copy.deepcopy(fieldFixed)
maskMovingCupy1 = copy.deepcopy(maskFixed)
maskMovingCupy3 = copy.deepcopy(maskFixed)
maskMovingSitk = copy.deepcopy(maskFixed)
print('Moving with translateData')
translateData(movingCupy1, translationInMM=translation, outputBox=outputBox, fillValue=-1000, tryGPU=True, interpOrder=1)
translateData(fieldMovingCupy1, translationInMM=translation, outputBox=outputBox, tryGPU=True, interpOrder=1)
translateData(maskMovingCupy1, translationInMM=translation, outputBox=outputBox, tryGPU=True, interpOrder=1)
translateData(movingCupy3, translationInMM=translation, outputBox=outputBox, fillValue=-1000, tryGPU=True, interpOrder=3)
translateData(fieldMovingCupy3, translationInMM=translation, outputBox=outputBox, tryGPU=True, interpOrder=3)
translateData(maskMovingCupy3, translationInMM=translation, outputBox=outputBox, tryGPU=True, interpOrder=3)
translateData(movingSitk, translationInMM=translation, outputBox=outputBox, fillValue=-1000)
translateData(fieldMovingSitk, translationInMM=translation, outputBox=outputBox)
translateData(maskMovingSitk, translationInMM=translation, outputBox=outputBox)
showImagesAndFieldAndMask(fixed, movingCupy1, movingCupy3, movingSitk, fieldFixed, fieldMovingCupy1,
fieldMovingCupy3, fieldMovingSitk,
maskFixed, maskMovingCupy1, maskMovingCupy3, maskMovingSitk, y_slice,
figTitle='Test using translateData', showField=showField, showImage=showImage,
showMask=showMask)
Test using rotateData#
print('-' * 40)
movingCupy1 = copy.deepcopy(fixed)
movingCupy3 = copy.deepcopy(fixed)
movingSitk = copy.deepcopy(fixed)
fieldMovingCupy1 = copy.deepcopy(fieldFixed)
fieldMovingCupy3 = copy.deepcopy(fieldFixed)
fieldMovingSitk = copy.deepcopy(fieldFixed)
maskMovingCupy1 = copy.deepcopy(maskFixed)
maskMovingCupy3 = copy.deepcopy(maskFixed)
maskMovingSitk = copy.deepcopy(maskFixed)
print('Moving with rotateData')
rotateData(movingCupy1, rotAnglesInDeg=rotation, outputBox=outputBox, fillValue=-1000, rotCenter=rotCenter, tryGPU=True, interpOrder=1)
rotateData(fieldMovingCupy1, rotAnglesInDeg=rotation, outputBox=outputBox, rotCenter=rotCenter, tryGPU=True, interpOrder=1)
rotateData(maskMovingCupy1, rotAnglesInDeg=rotation, outputBox=outputBox, rotCenter=rotCenter, tryGPU=True, interpOrder=1)
rotateData(movingCupy3, rotAnglesInDeg=rotation, outputBox=outputBox, fillValue=-1000, rotCenter=rotCenter, tryGPU=True, interpOrder=3)
rotateData(fieldMovingCupy3, rotAnglesInDeg=rotation, outputBox=outputBox, rotCenter=rotCenter, tryGPU=True, interpOrder=3)
rotateData(maskMovingCupy3, rotAnglesInDeg=rotation, outputBox=outputBox, rotCenter=rotCenter, tryGPU=True, interpOrder=3)
rotateData(movingSitk, rotAnglesInDeg=rotation, outputBox=outputBox, fillValue=-1000, rotCenter=rotCenter)
rotateData(fieldMovingSitk, rotAnglesInDeg=rotation, outputBox=outputBox, rotCenter=rotCenter)
rotateData(maskMovingSitk, rotAnglesInDeg=rotation, outputBox=outputBox, rotCenter=rotCenter)
showImagesAndFieldAndMask(fixed, movingCupy1, movingCupy3, movingSitk, fieldFixed, fieldMovingCupy1,
fieldMovingCupy3, fieldMovingSitk,
maskFixed, maskMovingCupy1, maskMovingCupy3, maskMovingSitk, y_slice,
figTitle='Test using rotateData', showField=showField, showImage=showImage,
showMask=showMask)