# Registration of CT Images with Large Organ Deformations

**URL:** <https://discourse.itk.org/t/registration-of-ct-images-with-large-organ-deformations/7627>\
**Category:** Beginner Questions\
**Tags:** dicom, simpleitk, registration, python\
**Created:** [September 2, 2025, 4:15pm UTC](https://discourse.itk.org/t/registration-of-ct-images-with-large-organ-deformations/7627 "2025-09-02T16:15:48Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![blue\_sky](https://discourse.itk.org/letter_avatar_proxy/v4/letter/b/dbc845/32.png) [@blue\_sky](https://discourse.itk.org/u/blue_sky)\
**Post date:** [September 2, 2025, 4:15pm UTC](https://discourse.itk.org/t/registration-of-ct-images-with-large-organ-deformations/7627/1 "2025-09-02T16:15:49Z")

</div>

Hello everyone,

I have two CT images that I want to register together. One is from external beam radiotherapy (EBRT), and the other is from brachytherapy. Before registration, I harmonized both images by making their Origin, dimensions, Spacing, and etc. identical. Additionally, I modified the Hounsfield units of the bladder and rectum in both images. Since the applicator is not present in the EBRT image, I set its Hounsfield unit to zero in the brachytherapy image.

After that, I ran the following code for registration. In writing parts of this code, I also drew from examples in the link: [Examples — SimpleITK 2.4.0 documentation](https://simpleitk.readthedocs.io/en/master/link_examples.html).  
However, despite the long execution time, I didn’t achieve good results. here is my result image. As you can see in the image, the registered image is not satisfactory.  
Can you guide me on how to register two CT images with significant changes in organs like the bladder and rectum, so that I can achieve good accuracy without the process taking too long?

Thanks a lot

Here is my code:

> import SimpleITK as sitk  
> import os  
> import time  
> import gc  
> sitk.ProcessObject.SetGlobalDefaultNumberOfThreads(0)  
> def read\_dicom\_series(directory):  
> print(f"Loading DICOM series from: {directory}“)  
> series\_IDs = sitk.ImageSeriesReader.GetGDCMSeriesIDs(directory)  
> if not series\_IDs:  
> raise ValueError(f"No DICOM series found in directory: {directory}”)  
> series\_files = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(directory, series\_IDs[0])  
> reader = sitk.ImageSeriesReader()  
> reader.SetFileNames(series\_files)  
> image = reader.Execute()  
> return sitk.Cast(image, sitk.sitkFloat32)  
> #Configuration  
> fixed\_ct\_dir = r’D:\fr1\CT\_Unified’  
> moving\_ct\_dir = r’D:\extr\CT\_Unified’  
> output\_dir = r’D:\registration\_output’  
> #Main Script  
> os.makedirs(output\_dir, exist\_ok=True)  
> dvf\_output\_path = os.path.join(output\_dir, ‘final\_dvf\_extr\_to\_fr1.mha’)  
> resampled\_moving\_ct\_path = os.path.join(output\_dir, ‘resampled\_moving\_ct.nii.gz’)  
> total\_start\_time = time.time()  
> fixed\_image = read\_dicom\_series(fixed\_ct\_dir)  
> moving\_image = read\_dicom\_series(moving\_ct\_dir)  
> print(f"Images loaded successfully. Fixed image size: {fixed\_image.GetSize()}“)  
> sampling\_percentage = 0.05  
> #Stages 1 & 2: Rigid and Affine  
> print(”\n— Starting Stage 1: Rigid Registration —“)  
> R\_rigid = sitk.ImageRegistrationMethod()  
> R\_rigid.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)  
> R\_rigid.SetMetricSamplingStrategy(R\_rigid.RANDOM)  
> R\_rigid.SetMetricSamplingPercentage(sampling\_percentage)  
> R\_rigid.SetInterpolator(sitk.sitkLinear)  
> R\_rigid.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=200, convergenceMinimumValue=1e-6, convergenceWindowSize=20)  
> R\_rigid.SetOptimizerScalesFromPhysicalShift()  
> R\_rigid.SetShrinkFactorsPerLevel(shrinkFactors=[8, 4, 2])  
> R\_rigid.SetSmoothingSigmasPerLevel(smoothingSigmas=[4, 2, 1])  
> R\_rigid.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()  
> initial\_rigid\_transform = sitk.CenteredTransformInitializer(fixed\_image, moving\_image, sitk.Euler3DTransform(), sitk.CenteredTransformInitializerFilter.GEOMETRY)  
> R\_rigid.SetInitialTransform(initial\_rigid\_transform)  
> rigid\_transform = R\_rigid.Execute(fixed\_image, moving\_image)  
> print(f"Rigid registration finished. Final Metric: {R\_rigid.GetMetricValue():.4f}”)  
> print(“\n— Starting Stage 2: Affine Registration —”)  
> R\_affine = sitk.ImageRegistrationMethod()  
> R\_affine.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)  
> R\_affine.SetMetricSamplingStrategy(R\_affine.RANDOM)  
> R\_affine.SetMetricSamplingPercentage(sampling\_percentage)  
> R\_affine.SetInterpolator(sitk.sitkLinear)  
> R\_affine.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=200, convergenceMinimumValue=1e-6, convergenceWindowSize=20)  
> R\_affine.SetOptimizerScalesFromPhysicalShift()  
> R\_affine.SetShrinkFactorsPerLevel(shrinkFactors=[4, 2])  
> R\_affine.SetSmoothingSigmasPerLevel(smoothingSigmas=[2, 1])  
> R\_affine.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()  
> initial\_affine\_transform = sitk.AffineTransform(3)  
> initial\_affine\_transform.SetMatrix(rigid\_transform.GetMatrix())  
> initial\_affine\_transform.SetTranslation(rigid\_transform.GetTranslation())  
> initial\_affine\_transform.SetCenter(rigid\_transform.GetCenter())  
> R\_affine.SetInitialTransform(initial\_affine\_transform)  
> affine\_transform = R\_affine.Execute(fixed\_image, moving\_image)  
> print(f"Affine registration finished. Final Metric: {R\_affine.GetMetricValue():.4f}“)  
> #Stage 3: Multi-level B-Spline Registration  
> print(”\n— Starting Stage 3: Multi-Level B-Spline Registration —“)  
> moving\_prealigned = sitk.Resample(moving\_image, fixed\_image, affine\_transform, sitk.sitkLinear, 0.0, moving\_image.GetPixelID())  
> #B-Spline Coarse Stage  
> print(”— Sub-stage: Coarse B-Spline —“)  
> R\_bspline\_coarse = sitk.ImageRegistrationMethod()  
> R\_bspline\_coarse.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)  
> R\_bspline\_coarse.SetMetricSamplingStrategy(R\_bspline\_coarse.RANDOM)  
> R\_bspline\_coarse.SetMetricSamplingPercentage(sampling\_percentage)  
> R\_bspline\_coarse.SetInterpolator(sitk.sitkLinear)  
> R\_bspline\_coarse.SetOptimizerAsLBFGSB(gradientConvergenceTolerance=1e-5, numberOfIterations=100, maximumNumberOfCorrections=5)  
> transform\_domain\_mesh\_size\_coarse = [10] \* fixed\_image.GetDimension()  
> R\_bspline\_coarse.SetInitialTransform(sitk.BSplineTransformInitializer(image1=fixed\_image, transformDomainMeshSize=transform\_domain\_mesh\_size\_coarse))  
> R\_bspline\_coarse.SetShrinkFactorsPerLevel([4, 2])  
> R\_bspline\_coarse.SetSmoothingSigmasPerLevel([2, 1])  
> bspline\_transform\_coarse = R\_bspline\_coarse.Execute(fixed\_image, moving\_prealigned)  
> print(f"B-Spline (Coarse) finished. Metric: {R\_bspline\_coarse.GetMetricValue():.4f}”)  
> #B-Spline Fine Stage  
> print(“— Sub-stage: Fine B-Spline —”)  
> R\_bspline\_fine = sitk.ImageRegistrationMethod()  
> R\_bspline\_fine.SetMetricAsANTSNeighborhoodCorrelation(radius=5)  
> R\_bspline\_fine.SetMetricSamplingStrategy(R\_bspline\_fine.RANDOM)  
> R\_bspline\_fine.SetMetricSamplingPercentage(sampling\_percentage)  
> R\_bspline\_fine.SetInterpolator(sitk.sitkLinear)  
> R\_bspline\_fine.SetOptimizerAsLBFGSB(gradientConvergenceTolerance=1e-5, numberOfIterations=100, maximumNumberOfCorrections=5)  
> transform\_domain\_mesh\_size\_fine = [12] \* fixed\_image.GetDimension()  
> R\_bspline\_fine.SetInitialTransform(bspline\_transform\_coarse, inPlace=False)  
> R\_bspline\_fine.SetShrinkFactorsPerLevel([2, 1])  
> R\_bspline\_fine.SetSmoothingSigmasPerLevel([1, 0])  
> bspline\_transform\_fine = R\_bspline\_fine.Execute(fixed\_image, moving\_prealigned)  
> print(f"B-Spline (Fine) finished. Final Metric: {R\_bspline\_fine.GetMetricValue():.4f}“)  
> #Final Stage: Combine and Save  
> print(”\n— Final Stage: Combining transforms and generating DVF —“)  
> final\_transform = sitk.CompositeTransform(affine\_transform)  
> final\_transform.AddTransform(bspline\_transform\_fine)  
> displacement\_filter = sitk.TransformToDisplacementFieldFilter()  
> displacement\_filter.SetReferenceImage(fixed\_image)  
> final\_dvf = displacement\_filter.Execute(final\_transform)  
> final\_dvf = sitk.Cast(final\_dvf, sitk.sitkVectorFloat64)  
> print(“Clearing intermediate objects from memory…”)  
> del R\_rigid, R\_affine, R\_bspline\_coarse, R\_bspline\_fine, rigid\_transform, affine\_transform, bspline\_transform\_coarse, bspline\_transform\_fine, moving\_prealigned, final\_transform, displacement\_filter  
> gc.collect()  
> print(“Writing final DVF…”)  
> sitk.WriteImage(final\_dvf, dvf\_output\_path)  
> print(f"Final DVF saved to: {dvf\_output\_path}”)  
> print(“\n— Generating resampled CT —”)  
> dvf\_transform = sitk.DisplacementFieldTransform(final\_dvf)  
> resampled\_moving = sitk.Resample(moving\_image, fixed\_image, dvf\_transform, sitk.sitkLinear, 0.0)  
> sitk.WriteImage(resampled\_moving, resampled\_moving\_ct\_path)  
> print(f"Final resampled moving CT (for validation) saved to: {resampled\_moving\_ct\_path}“)  
> print(f”\nRegistration completed successfully in {time.time() - total\_start\_time:.2f} seconds.")

 ![registration result](https://discourse.itk.org/uploads/default/original/2X/5/5228b50b62d3a454f9bec8144219f55c30a1991d.jpeg)

---

<div class="post-metadata">

**Author:** ![zivy](https://discourse.itk.org/user_avatar/discourse.itk.org/zivy/32/1726_2.png) [@zivy](https://discourse.itk.org/u/zivy)\
**Post date:** [September 2, 2025, 7:49pm UTC](https://discourse.itk.org/t/registration-of-ct-images-with-large-organ-deformations/7627/2 "2025-09-02T19:49:17Z")

</div>

Hello @blue_sky,

The changes made to the images prior to registration, origin/spacing/etc., are not necessary. You should leave the original scans as is.

This “harmonization” step is introducing unnecessary resampling which may have a detrimental effect on registration. The modifications of the intensities for bladder and rectum are also likely not helpful. Finally, setting the applicator pixels to zero is likely not what you want. A HU value of 0 corresponds to water, which is likely the same value as urine in the bladder. If you want to change the intensity of the applicator a value of -1000, air, is probably more appropriate, but I would leave it as is.

Before trying to do this programmatically, please try using Slicer (see [registration documentation](https://slicer.readthedocs.io/en/latest/user_guide/registration.html#automatic-image-registration)). The interface for the advanced registration is friendlier than the programmatic SimpleITK interface and the default task specific values for the optimization settings mean you do not need expertise in registration to set them in an appropriate way for your setting (CT/CT non-rigid registration).