All the examples might not be up to date on the website, please go to our GitLab repository and download latest version of the notebook.

Rigid registration

import copy
import os
import numpy as np
import matplotlib.pyplot as plt
import time
import logging

from opentps.core.processing.registration.registrationRigid import RegistrationRigid
from opentps.core.examples.syntheticData import *
from opentps.core.processing.imageProcessing.resampler3D import resampleImage3DOnImage3D
from opentps.core.processing.imageProcessing.imageTransform3D import rotateData, translateData

Generate synthetic images

fixed = createSynthetic3DCT()
moving = copy.copy(fixed)

translation = np.array([15, 0, 10])
rotation = np.array([0, 5, 0])

translateData(moving, translation, outputBox='same')
rotateData(moving, rotation, outputBox='same')

Perform the Rigid registration

start_time = time.time()
reg = RegistrationRigid(fixed, moving)
transform = reg.compute()

processing_time = time.time() - start_time
print('Registration processing time was', processing_time, '\n')
print('Translation', transform.getTranslation())
print('Rotation in deg', transform.getRotationAngles(inDegrees=True), '\n')

## Two ways of getting the deformed moving image
#deformedImage = reg.deformed
deformedImage = transform.deformImage(moving)

## Resample it to the same grid as the fixed image
resampledOnFixedGrid = resampleImage3DOnImage3D(deformedImage, fixedImage=fixed, fillValue=-1000)

Final metric value: 19945.711049659527
Optimizer's stopping condition, RegularStepGradientDescentOptimizerv4: Step too small after 31 iterations. Current step (9.53674e-07) is less than minimum step (1e-06).
itk::simple::CompositeTransform
 CompositeTransform (0000013E3A2C9380)
   RTTI typeinfo:   class itk::CompositeTransform<double,3>
   Reference Count: 1
   Modified Time: 1330522
   Debug: Off
   Object Name: 
   Observers: 
     none
   Transforms in queue, from begin to end:
   >>>>>>>>>
   Euler3DTransform (0000013E1D5AF610)
     RTTI typeinfo:   class itk::Euler3DTransform<double>
     Reference Count: 1
     Modified Time: 1330512
     Debug: Off
     Object Name: 
     Observers: 
       none
     Matrix: 
       0.995794 -0.000422025 0.0916217 
       0.000452553 1 -0.000312418 
       -0.0916215 0.000352568 0.995794 
     Offset: [14.9122, -0.000670789, 1.13981]
     Center: [84.5, 84.5, 99]
     Translation: [23.5916, 0.00662771, -6.98882]
     Inverse: 
       0.995794 0.000452553 -0.0916215 
       -0.000422025 1 0.000352568 
       0.0916217 -0.000312418 0.995794 
     Singular: 0
     Euler's angles: AngleX=0.000352568 AngleY=0.0917502 AngleZ=0.000422025
     m_ComputeZYX = 0
   End of MultiTransform.
<<<<<<<<<<
   TransformsToOptimizeFlags, begin() to end(): 
      1 
   TransformsToOptimize in queue, from begin to end:
   End of TransformsToOptimizeQueue.
<<<<<<<<<<
   End of CompositeTransform.
<<<<<<<<<<

Registration processing time was 3.0145018100738525 

Translation [-2.35916332e+01 -6.62771221e-03  6.98882239e+00]
Rotation in deg [-0.02020065 -5.25689914 -0.02418025] 

Compute the difference between the 2 images and check results of the registration

#IMAGE DIFFERENCE
diff_before = fixed.copy()
diff_before._imageArray = fixed.imageArray - moving.imageArray
diff_after = fixed.copy()
diff_after._imageArray = fixed.imageArray - resampledOnFixedGrid.imageArray

# CHECK RESULTS
diff_before_sum = abs(diff_before.imageArray).sum()
diff_after_sum = abs(diff_after.imageArray).sum()
assert diff_before_sum - diff_after_sum > 0, f"Image difference is larger after registration"

Plot results

#Output path
output_path = 'Output'
if not os.path.exists(output_path):
    os.makedirs(output_path)

y_slice = 95
fig, ax = plt.subplots(2, 3,figsize=(8,10))
ax[0, 0].imshow(fixed.imageArray[:, y_slice, :],cmap='gray')
ax[0, 0].set_title('Fixed')
ax[0, 0].set_xlabel('Origin: '+f'{fixed.origin[0]}'+','+f'{fixed.origin[1]}'+','+f'{fixed.origin[2]}')
ax[0, 1].imshow(moving.imageArray[:, y_slice, :],cmap='gray')
ax[0, 1].set_title('Moving')
ax[0, 1].set_xlabel('Origin: ' + f'{moving.origin[0]}' + ',' + f'{moving.origin[1]}' + ',' + f'{moving.origin[2]}')
diffBef = ax[0, 2].imshow(abs(diff_before.imageArray[:, y_slice, :]), vmin=0, vmax=1500)
ax[0, 2].set_title('Diff before')
fig.colorbar(diffBef, ax=ax[0, 2])
ax[1, 0].imshow(deformedImage.imageArray[:, y_slice, :],cmap='gray')
ax[1, 0].set_title('DeformedMoving')
ax[1, 0].set_xlabel('Origin: ' + f'{deformedImage.origin[0]:.1f}' + ',' + f'{deformedImage.origin[1]:.1f}' + ',' + f'{deformedImage.origin[2]:.1f}')
ax[1, 1].imshow(resampledOnFixedGrid.imageArray[:, y_slice, :],cmap='gray')
ax[1, 1].set_title('resampledOnFixedGrid')
ax[1, 1].set_xlabel('Origin: ' + f'{resampledOnFixedGrid.origin[0]:.1f}' + ',' + f'{resampledOnFixedGrid.origin[1]:.1f}' + ',' + f'{resampledOnFixedGrid.origin[2]:.1f}')
diffAft = ax[1, 2].imshow(abs(diff_after.imageArray[:, y_slice, :]), vmin=0, vmax=1500)
ax[1, 2].set_title('Diff after')
#fig.colorbar(diffAft, ax=ax[1, 2])
plt.savefig(os.path.join(output_path,"Rigid1.png"),format = 'png')
plt.close()

png

Download this notebook via our GitLab repository