VarianAPIs / VelocityEngine

Python and C# Clients for Velocity Scripting API
7 stars 5 forks source link

Is there a random component in deformable registration? #21

Closed mfizyczka closed 2 years ago

mfizyczka commented 2 years ago

I used this code to generate 5 times the same registration and propagate contours. Results are different with each run. Why?

for i in range(1,6):
  print(i)
  deformableRegistrationSettings = velocity.BSplineDeformableRegistrationSettingsStructure()
  # create registration and perform registration:
  registrationName = 'apiDefMP' + str(i)
  registration = regOps.createNewRegistration(registrationName)
  orThrow(registration, regOps)
  orThrow(engine.loadRegistration(registration.getVelocityId()))
  orThrow(regOps.performBsplineRegistrationDICOM(deformableRegistrationSettings), regOps)
  registration = regOps.saveRegistration()
  orThrow(registration, regOps)
  # create a new structure set on the primary volume:
  targetSetName = registrationName
  targetSet = strOps.createStructureSet(targetSetName, True)
  orThrow(targetSet, strOps)
  # propagate structures:
  orThrow(strOps.copyStructuresToPrimary([ctv5425.getVelocityId(), ctv7000.getVelocityId()], targetSet.getVelocityId()), strOps)
  # save set:
  orThrow(strOps.saveStructureSet(targetSet.getVelocityId()), strOps)

Axial113_random DefRegistrationMP

Performing rigid registration first (with auto alignment) gives consistent results: Axial113_AutoAlignRigidMP_api

However in this case I had to run code 5 times. With loop I've got an error. I don't know how to correct my code.

i = 5
registrationName = 'apiDfSetRMP' + str(i)
registration = regOps.createNewRegistration(registrationName)
# create registration and perform registration:
orThrow(registration, regOps)
orThrow(engine.loadRegistration(registration.getVelocityId()))
print('Running rigid registration...')
rigidDefaultRegistrationSettings = velocity.DefaultRigidRegistrationSettings()
orThrow(regOps.performRigidRegistrationDICOM(rigidDefaultRegistrationSettings), regOps)
print('Performing deformable registration...')
deformableDefaultRegistrationSettings = velocity.DefaultBSplineDeformableRegistrationSettings()
orThrow(regOps.performBsplineRegistrationDICOM(deformableDefaultRegistrationSettings), regOps)
registration = regOps.saveRegistration()
orThrow(registration, regOps)
# create a new structure set on the primary volume:
targetSetName = registrationName
targetSet = strOps.createStructureSet(targetSetName, True)
orThrow(targetSet, strOps)
# propagate structures:
orThrow(strOps.copyStructuresToPrimary([ctv5425.getVelocityId(), ctv7000.getVelocityId()], targetSet.getVelocityId()), strOps)
# save set:
orThrow(strOps.saveStructureSet(targetSet.getVelocityId()), strOps)

Output:

Running rigid registration...
Performing deformable registration...
Running rigid registration...
Traceback (most recent call last):
  File "d:\IMAGINATION\#VelocityEngine\PythonScripting\Registrations\testAndLearnAboutDeformableRegistrationSettings.py", line 290, in <module>
    orThrow(regOps.performRigidRegistrationDICOM(rigidRegistrationSettings), regOps)
  File "d:\IMAGINATION\#VelocityEngine\PythonScripting\Registrations\testAndLearnAboutDeformableRegistrationSettings.py", line 36, in orThrow
    raise RuntimeError(e.getErrorMessage())
RuntimeError: Unable to perform rigid registration because a deformable object is currently set
jakecobb commented 2 years ago

As with the Velocity GUI, when you create a new registration it starts as a copy of the currently loaded registration. So when you run it repeatedly as shown here you are not starting from the same state, but instead each run is starting with the previous result.

Internally there is always a registration of some sort loaded. In the UI, if no registration is checked, the internal registration is such that the centers of the two volumes are aligned. The API is currently missing a function to let you revert to this "no registration", you can either load the one you want to copy prior to createRegistration or unloading and reloading the volumes will also clear the registration.

Running deformation on a rigid registration converts it to a deformable (which includes a rigid component), but you can't run rigid on a deformable to convert back to a rigid.

afwaller commented 2 years ago

Note: image registration in Velocity is deterministic and will yield the same result on the same system every time.

There are minor changes to improve registration (aspects such as performance, memory usage, bug fixes, noise handling, multi threading) between versions and releases.

Additionally, depending on the system and processor, slightly different results can occur between various computer systems with different processors and processor counts, depending on how the registration is actually performed. That said, image registration is, as previously noted, deterministic for a given version and system.

mfizyczka commented 2 years ago

I still see somehow different results with API:

def loadVolumesAndStructures(pCtName, rCtName):
  # load primary volume - volume with originally delineated structures
  primaryVolumeUID = patientVolumes[rCtName][1]
  orThrow(engine.loadPrimaryVolumeByUID(primaryVolumeUID))
  global primaryVolume
  primaryVolume = engine.getPrimaryVolume()

  # load secondary volume - to this volume structures will be progressed
  secondaryVolumeUID = patientVolumes[pCtName][1]
  orThrow(engine.loadSecondaryVolumeByUID(secondaryVolumeUID))
  global secondaryVolume
  secondaryVolume = engine.getSecondaryVolume()

  # create dictionary of StructureSets on primaryVolume
  secondaryVolumeStructureSets = {}
  for ss in secondaryVolume.getStructureSets(): secondaryVolumeStructureSets[ss.getName()] = [ss.getVelocityId(), ss.getInstanceUID()]

  # load secondary set and structures:
  secondarySetName = 'RS: Unapproved'
  secondarySet  = next(s for s in secondaryVolume.getStructureSets() if s.getName() == secondarySetName)
  global ctv7000
  ctv7000 = next(s for s in secondarySet.getStructures() if s.getName() == 'CTV_7000')
  global ctv5425
  ctv5425 = next(s for s in secondarySet.getStructures() if s.getName() == 'CTV_5425')

  return True

 DEFAULT REGISTRATION SETTINGS (Multi Pass)
for i in range(1,4):
  # load Volumes:
  loadVolumesAndStructures('CT[Accession N.52]', 'CT[Accession N.58]')
  deformableDefaultRegistrationSettings = velocity.DefaultBSplineDeformableRegistrationSettings()
  # create registration and perform registration:
  registrationName = 'apiDefSetMP_' + str(i)
  print(registrationName)
  registration = regOps.createNewRegistration(registrationName)
  orThrow(registration, regOps)
  orThrow(engine.loadRegistration(registration.getVelocityId()))
  orThrow(regOps.performBsplineRegistrationDICOM(deformableDefaultRegistrationSettings), regOps)
  registration = regOps.saveRegistration()
  orThrow(registration, regOps)
  # create a new structure set on the primary volume:
  targetSetName = registrationName
  targetSet = strOps.createStructureSet(targetSetName, True)
  orThrow(targetSet, strOps)
  # propagate structures:
  orThrow(strOps.copyStructuresToPrimary([ctv5425.getVelocityId(), ctv7000.getVelocityId()], targetSet.getVelocityId()), strOps)
  # save set:
  orThrow(strOps.saveStructureSet(targetSet.getVelocityId()), strOps)
  # unload volumes:
  orThrow(engine.unloadPrimaryVolume, engine)
  orThrow(engine.unloadSecondaryVolume, engine)

Comparison in GUI: apiDefSetMP_1_2_3

jakecobb commented 2 years ago

You have not unloaded the volumes because you passed the functions instead of calling them:

  # unload volumes:
  orThrow(engine.unloadPrimaryVolume(), engine)
  orThrow(engine.unloadSecondaryVolume(), engine)
mfizyczka commented 2 years ago

Oh... It's true. My mistake.