# Unifying CT Images Geometry and Cropping to Common ROI
**URL:** https://discourse.itk.org/t/unifying-ct-images-geometry-and-cropping-to-common-roi/7614
**Category:** Beginner Questions
**Tags:** registration, python, github, simpleitk
**Created:** [August 15, 2025, 10:59am UTC](https://discourse.itk.org/t/unifying-ct-images-geometry-and-cropping-to-common-roi/7614 "2025-08-15T10:59:38Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@blue\_sky](https://discourse.itk.org/u/blue_sky)
#### Post date: [August 15, 2025, 10:59am UTC](https://discourse.itk.org/t/unifying-ct-images-geometry-and-cropping-to-common-roi/7614/1 "2025-08-15T10:59:38Z")
Hi everyone,
I am currently working on a medical image processing task that involves the registration of five CT image series for a single patient. My dataset consists of:
One CT scan acquired during an External Radiotherapy (EBRT) session, which presents a significantly different Field of View compared to the others and four subsequent CT scans from Brachytherapy for the cervix.
To continue my work, a crucial pre-processing step is required:
- Homogenizing the spatial information (Origin, Image Spacing, and FOV) across all five images to ensure they reside in a consistent physical space.
- Extracting the largest common anatomical region that is present in all scans, to serve as the input for subsequent registration and analysis steps.
I have been exploring the SimpleITK-Notebooks for guidance, particularly the “Create reference domain” section in:
([SimpleITK-Notebooks/Python/70\_Data\_Augmentation.ipynb at main · InsightSoftwareConsortium/SimpleITK-Notebooks · GitHub](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/main/Python/70_Data_Augmentation.ipynb)) to define a unified image space.
Despite this, I am consistently encountering persistent and undesirable padding in various anatomical views (axial, sagittal, and coronal) of the processed images. Furthermore, some of my attempts have resulted in excessive cropping, inadvertently removing valuable anatomical structures required for downstream analysis.
Could you please provide guidance on how to robustly achieve a clean, common anatomical region across all these diverse CT scans, ensuring no residual padding remains while preserving essential anatomical features (like applicators or the bladder) and avoiding over-cropping?
Thanks a lot
my code:
> import SimpleITK as sitk
> import numpy as np
> import os
> def read\_dicom\_series(directory):
> “”"
> Reads a DICOM series from a directory and returns a SimpleITK Image object.
> “”"
> if not os.path.isdir(directory):
> raise FileNotFoundError(f"Directory not found: {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()
> print(f"Loaded DICOM series from {directory} with size {image.GetSize()}, origin {image.GetOrigin()}, spacing {image.GetSpacing()}“)
> return image
> def compute\_physical\_bounds(img):
> “””
> Computes the physical bounds (min and max coordinates) of an image in physical space.
> Returns a tuple of (min\_bounds, max\_bounds) for each dimension.
> “”"
> size = np.array(img.GetSize())
> spacing = np.array(img.GetSpacing())
> origin = np.array(img.GetOrigin())
> max\_bounds = origin + (size - 1) \* spacing
> return origin, max\_bounds
> def crop\_to\_common\_region(images):
> “”"
> Crops all images to the common anatomical region (intersection of physical bounds).
> Returns a list of cropped images.
> “”"
> if not images:
> raise ValueError(“No images provided for cropping.”)
> # Computes physical bounds for all images
> bounds = [compute\_physical\_bounds(img) for img in images]
> min\_bounds\_all = [b[0] for b in bounds] # List of origins
> max\_bounds\_all = [b[1] for b in bounds] # List of max coordinates
> # FindS the common region (max of min bounds, min of max bounds)
> common\_min = np.max(min\_bounds\_all, axis=0)
> common\_max = np.min(max\_bounds\_all, axis=0)
> if np.any(common\_min \>= common\_max):
> print(“Warning: No common region found among the images. Returning original images.”)
> return images
> print(f"Common physical region for cropping: Min bounds {common\_min}, Max bounds {common\_max}“)
> # Crops each image to the common region
> cropped\_images = []
> for i, img in enumerate(images):
> origin = np.array(img.GetOrigin())
> spacing = np.array(img.GetSpacing())
> size = np.array(img.GetSize())
> # Compute index bounds for cropping
> start\_idx = np.maximum(0, np.ceil((common\_min - origin) / spacing)).astype(int)
> end\_idx = np.minimum(size, np.floor((common\_max - origin) / spacing) + 1).astype(int)
> if np.any(end\_idx \<= start\_idx):
> print(f"Warning: Cannot crop image {i+1} (invalid bounds {start\_idx} to {end\_idx}). Keeping original.”)
> cropped\_images.append(img)
> continue
> # Extract region
> cropped\_img = img[start\_idx[0]:end\_idx[0], start\_idx[1]:end\_idx[1], start\_idx[2]:end\_idx[2]]
> print(f"Cropped image {i+1} to size {cropped\_img.GetSize()}“)
> cropped\_images.append(cropped\_img)
> return cropped\_images
> def create\_reference\_domain(images, pixel\_size\_per\_dimension=128):
> “””
> Uses the largest physical size among images.
> “”"
> dimension = images[0].GetDimension()
> # Computes the largest physical size among all images
> reference\_physical\_size = np.zeros(dimension)
> for img in images:
> phys\_size = np.array([(sz - 1) \* spc for sz, spc in zip(img.GetSize(), img.GetSpacing())])
> reference\_physical\_size = np.maximum(reference\_physical\_size, phys\_size)
> print(f"Reference physical size (largest FOV): {reference\_physical\_size}“)
> # Sets origin and direction to standard values (as in the notebook)
> reference\_origin = np.zeros(dimension)
> reference\_direction = np.identity(dimension).flatten()
> # Sets arbitrary pixel size
> reference\_size = [pixel\_size\_per\_dimension] \* dimension
> reference\_spacing = [phys\_sz / (sz - 1) for sz, phys\_sz in zip(reference\_size, reference\_physical\_size)]
> print(f"Reference size: {reference\_size}”)
> print(f"Reference spacing: {reference\_spacing}“)
> # Creates reference image
> reference\_image = sitk.Image(reference\_size, images[0].GetPixelIDValue())
> reference\_image.SetOrigin(reference\_origin)
> reference\_image.SetSpacing(reference\_spacing)
> reference\_image.SetDirection(reference\_direction)
> # Computes reference center (as in the notebook)
> reference\_center = np.array(
> reference\_image.TransformContinuousIndexToPhysicalPoint(
> np.array(reference\_image.GetSize()) / 2.0
> )
> )
> print(f"Reference center: {reference\_center}”)
> return reference\_image
> def resample\_to\_reference(images, reference\_image):
> “”"
> Resamples all images to the reference domain, centering them.
> “”"
> resampled\_images = []
> for i, img in enumerate(images):
> # Creates a transform to center the image in the reference domain
> transform = sitk.CenteredTransformInitializer(
> reference\_image,
> img,
> sitk.Euler3DTransform(),
> sitk.CenteredTransformInitializerFilter.GEOMETRY
> )
> resampler = sitk.ResampleImageFilter()
> resampler.SetReferenceImage(reference\_image)
> resampler.SetInterpolator(sitk.sitkBSpline)
> resampler.SetTransform(transform)
> resampler.SetDefaultPixelValue(0)
> resampled\_img = resampler.Execute(img)
> print(f"Resampled image {i+1} to reference domain with size {resampled\_img.GetSize()}“)
> resampled\_images.append(resampled\_img)
> return resampled\_images
> dicom\_dirs = [
> r’C:\Users\Desktop\extrenal’,
> r’C:\Users\Desktop\fr1’,
> r’C:\Users\Desktop\fr2’,
> r’C:\Users\Desktop\fr3’,
> r’C:\Users\Desktop\fr4’
> ]
> output\_dir = r’C:\Users\Desktop\resampled’
> if not os.path.exists(output\_dir):
> os.makedirs(output\_dir)
> # Step 1: Loads DICOM series
> images = [read\_dicom\_series(dir\_path) for dir\_path in dicom\_dirs]
> # Step 2: Crops to common region (to preserve only the intersection of physical regions)
> do\_crop\_to\_common = True
> if do\_crop\_to\_common:
> processed\_images = crop\_to\_common\_region(images)
> else:
> processed\_images = images
> print(“Skipping cropping to common region; using original images directly.”)
> # Step 3: Creates reference domain (using largest physical size)
> reference\_image = create\_reference\_domain(processed\_images, pixel\_size\_per\_dimension=128)
> # Step 4: Resamples all images to reference domain with centering
> resampled\_images = resample\_to\_reference(processed\_images, reference\_image)
> # Step 5: Saves resampled images (as NIfTI)
> for i, resampled\_img in enumerate(resampled\_images):
> output\_path = os.path.join(output\_dir, f’resampled\_ct{i+1}.nii’)
> sitk.WriteImage(resampled\_img, output\_path)
> print(f"Saved resampled image to {output\_path}”)
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [August 15, 2025, 2:42pm UTC](https://discourse.itk.org/t/unifying-ct-images-geometry-and-cropping-to-common-roi/7614/2 "2025-08-15T14:42:03Z")
Hello @blue_sky,
Based on your goal, registering multiple subsequent CT scans to an initial scan, the data modification step that you are attempting to perform is best avoided. It introduces a resampling operation before any registration that modifies the original data unnecessarily.
Also, based on your description “One CT scan acquired during an External Radiotherapy (EBRT) session, which presents a significantly different Field of View compared to the others” the physical size of the reference domain returned by `create_reference_domain` is potentially very large and thus the resulting resampled images likely have a significant amount of empty space, which is what you are seeing. In the extreme, the first scan is a whole body scan and the consecutive scans are of the pelvic region. After this processing the pelvic region scans are embedded in a much larger image which is mostly empty.
Given the problem description, a more relevant workflow is the one illustrated in the [x-ray panorma notebook](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/main/Python/69_x-ray-panorama.ipynb). Please take a look at that flow, registration followed by computing the shared physical domain and resampling onto that. The relevant function which is called after registration is `create_images_in_shared_coordinate_system`.
If this does not address your needs, please provide additional details.
---
### Author:  [@blue\_sky](https://discourse.itk.org/u/blue_sky)
#### Post date: [August 22, 2025, 7:16am UTC](https://discourse.itk.org/t/unifying-ct-images-geometry-and-cropping-to-common-roi/7614/3 "2025-08-22T07:16:28Z")
hello, Many thanks for your guidance. I have implemented a workflow that:
-Rigidly registers CT images.
-Resamples the corresponding RTDOSE images to match the registered CT geometry.
-Modifies RTSTRUCT contours according to the CT transformation, ensuring they match perfectly with the updated CT.
Currently, the code saves the output CT, Dose, and Mask images in NIfTI format. When I load these results into 3D Slicer, they display correctly — spacing, origin, geometry, and registration are all preserved.
However, for the next steps in my workflow, I need to save these outputs in standard DICOM format:
CT images should be stored as a proper DICOM series.
Dose images should preserve geometry, DoseGridScaling, and match exactly to the CT’s spacing, origin, and size.
RTSTRUCT should remain spatially aligned with the transformed CT.
When I try saving as DICOM directly from the resampled SimpleITK images, I encounter issues:
Dose files don’t always preserve metadata correctly.
In some cases, resulting DICOMs are not readable or appear distorted in 3D Slicer.
Could you guide me on the best practices for converting these registered NIfTI outputs into fully compliant DICOM CT, DOSE, and RTSTRUCT files while preserving all geometry and metadata?
thanks a lot
My Code:
> import SimpleITK as sitk
> import numpy as np
> import os
> import pydicom
> from pydicom.uid import generate\_uid
> import logging
> from datetime import datetime
> import time
> def setup\_logging():
> “”“Setup logging configuration”“”
> logging.basicConfig(
> level=logging.INFO,
> format=‘%(asctime)s - %(levelname)s - %(message)s’,
> handlers=[
> logging.FileHandler(‘registration.log’),
> logging.StreamHandler()
> ]
> )
> return logging.getLogger( **name** )
> def read\_dicom\_series(directory, is\_dose=False):
> “”“Reads DICOM series with proper dose scaling”“”
> if not os.path.isdir(directory):
> raise FileNotFoundError(f"Directory not found: {directory}“)
> if is\_dose:
> dicom\_files = [os.path.join(directory, f) for f in os.listdir(directory)
> if os.path.isfile(os.path.join(directory, f)) and
> not f.lower().endswith((‘.txt’, ‘.xml’, ‘.log’, ‘.zip’))]
> if len(dicom\_files) != 1:
> raise ValueError(f"Expected ONE DICOM file in dose directory, but found {len(dicom\_files)} in: {directory}”)
> reader = sitk.ImageFileReader()
> reader.SetFileName(dicom\_files[0])
> image = reader.Execute()
> # Gets dose scaling factor
> dose\_scaling = 1.0
> try:
> dose\_scaling = float(reader.GetMetaData(“3004|000e”))
> print(f"Found and applied DoseGridScaling: {dose\_scaling}“)
> except (RuntimeError, ValueError):
> print(“Warning: Could not read DoseGridScaling tag. Assuming 1.0.”)
> return sitk.Cast(image, sitk.sitkFloat64) \* dose\_scaling, dicom\_files[0]
> else:
> series\_IDs = sitk.ImageSeriesReader.GetGDCMSeriesIDs(directory)
> if not series\_IDs:
> raise ValueError(f"No DICOM series found in: {directory}”)
> series\_files = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(directory, series\_IDs[0])
> reader = sitk.ImageSeriesReader()
> reader.SetFileNames(series\_files)
> image = reader.Execute()
> return image, series\_files[0]
> def select\_reference\_image(images, image\_names):
> “”“Selects reference image with smallest Z FOV”“”
> fovs\_z = [img.GetSize()[2] \* img.GetSpacing()[2] for img in images]
> min\_z\_fov\_index = np.argmin(fovs\_z)
> ref\_image = images[min\_z\_fov\_index]
> print(f"\nReference image selected: ‘{image\_names[min\_z\_fov\_index]}’ (Image {min\_z\_fov\_index + 1})“)
> return ref\_image, min\_z\_fov\_index
> def robust\_rigid\_registration(fixed\_image, moving\_image):
> “”“rigid registration””"
> # Converts to float32 for registration
> fixed\_image\_float = sitk.Cast(fixed\_image, sitk.sitkFloat32)
> moving\_image\_float = sitk.Cast(moving\_image, sitk.sitkFloat32)
> # Initialize transform
> initial\_transform = sitk.CenteredTransformInitializer(
> fixed\_image\_float, moving\_image\_float,
> sitk.Euler3DTransform(),
> sitk.CenteredTransformInitializerFilter.GEOMETRY
> )
> # Setups registration method
> R = sitk.ImageRegistrationMethod()
> R.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
> R.SetMetricSamplingStrategy(R.RANDOM)
> R.SetMetricSamplingPercentage(0.01)
> R.SetInterpolator(sitk.sitkLinear)
> R.SetOptimizerAsGradientDescent(
> learningRate=1.0,
> numberOfIterations=200,
> convergenceMinimumValue=1e-6,
> convergenceWindowSize=10
> )
> R.SetOptimizerScalesFromPhysicalShift()
> R.SetShrinkFactorsPerLevel(shrinkFactors=[4, 2, 1])
> R.SetSmoothingSigmasPerLevel(smoothingSigmas=[2, 1, 0])
> R.SetSmoothingSigmasAreSpecifiedInPhysicalUnits(True)
> R.SetInitialTransform(initial\_transform, inPlace=False)
> try:
> final\_transform = R.Execute(fixed\_image\_float, moving\_image\_float)
> print(f" Registration successful. Final metric: {R.GetMetricValue():.6f}“)
> return final\_transform
> except Exception as e:
> print(f” Registration failed: {e}. Returning initial alignment.“)
> return initial\_transform
> def create\_body\_mask(image, threshold=-500):
> “”“Creates body mask””"
> mask = sitk.BinaryThreshold(image, lowerThreshold=threshold, upperThreshold=10000,
> insideValue=1, outsideValue=0)
> return sitk.Cast(mask, sitk.sitkUInt8)
> def transform\_rt\_structure(rt\_file\_path, transform, reference\_image, output\_dir, dataset\_name):
> “”“Transforms RT Structure”“”
> try:
> # Reads RT Structure
> ds = pydicom.dcmread(rt\_file\_path)
> if not hasattr(ds, ‘ROIContourSequence’):
> logging.warning(f"No ROI contours found in {os.path.basename(rt\_file\_path)}“)
> return False
> # Gets inverse transform for structure transformation
> inverse\_transform = transform.GetInverse()
> roi\_count = 0
> contour\_count = 0
> # Transforms each ROI contour
> for roi\_contour in ds.ROIContourSequence:
> if hasattr(roi\_contour, ‘ContourSequence’):
> roi\_count += 1
> for contour in roi\_contour.ContourSequence:
> if hasattr(contour, ‘ContourData’):
> points = np.array(contour.ContourData).reshape(-1, 3)
> # Transforms points using inverse transform
> transformed\_points = []
> for point in points:
> transformed\_point = inverse\_transform.TransformPoint(point.tolist())
> transformed\_points.extend(transformed\_point)
> contour.ContourData = transformed\_points
> contour\_count += 1
> # Updates Frame of Reference UID to match reference image
> try:
> if hasattr(ds, ‘ReferencedFrameOfReferenceSequence’) and ds.ReferencedFrameOfReferenceSequence:
> ds.ReferencedFrameOfReferenceSequence[0].FrameOfReferenceUID = generate\_uid()
> except:
> pass
> # Generates new SOP Instance UID
> ds.SOPInstanceUID = generate\_uid()
> # Saves transformed RT Structure
> output\_path = os.path.join(output\_dir, f”{dataset\_name}\_RT\_Structure.dcm")
> ds.save\_as(output\_path)
> print(f" Successfully transformed RT Structure: {roi\_count} ROIs, {contour\_count} contours")
> print(f" Saved to: {os.path.basename(output\_path)}“)
> return True
> except Exception as e:
> print(f” ERROR transforming RT Structure {os.path.basename(rt\_file\_path)}: {str(e)}“)
> return False
> def main():
> “”“Main processing function””"
> logger = setup\_logging()
> # Directory paths
> base\_dir = r’C:\Users\resampled’
> ct\_dicom\_dirs = [
> os.path.join(base\_dir, r’ct\extr’),
> os.path.join(base\_dir, r’ct\fr1’),
> os.path.join(base\_dir, r’ct\fr2’),
> os.path.join(base\_dir, r’ct\fr3’),
> os.path.join(base\_dir, r’ct\fr4’)
> ]
> dose\_dicom\_dirs = [
> os.path.join(base\_dir, r’dose\extr’),
> os.path.join(base\_dir, r’dose\fr1’),
> os.path.join(base\_dir, r’doses\fr2’),
> os.path.join(base\_dir, r’doses\fr3’),
> os.path.join(base\_dir, r’doses\fr4’)
> ]
> rt\_struct\_dirs = [
> os.path.join(base\_dir, r’rts\extr’),
> os.path.join(base\_dir, r’rts\fr1’),
> os.path.join(base\_dir, r’rts\fr2’),
> os.path.join(base\_dir, r’rts\fr3’),
> os.path.join(base\_dir, r’rts\fr4’)
> ]
> output\_dir = r’C:\Users\resampled\_final’
> os.makedirs(output\_dir, exist\_ok=True)
> # Dataset names
> ct\_names = [‘extr’, ‘fr1’, ‘fr2’, ‘fr3’, ‘fr4’]
> print(“=== COMBINED REGISTRATION WORKFLOW ===”)
> print(“— Loading All Data —”)
> # Loads CT images
> ct\_images = []
> ct\_ref\_files = []
> for i, ct\_dir in enumerate(ct\_dicom\_dirs):
> print(f"Loading CT from: {os.path.basename(ct\_dir)}“)
> ct\_image, ref\_file = read\_dicom\_series(ct\_dir)
> ct\_images.append(ct\_image)
> ct\_ref\_files.append(ref\_file)
> print(f” CT size: {ct\_image.GetSize()}“)
> # Loads dose images
> dose\_images = []
> dose\_ref\_files = []
> for i, dose\_dir in enumerate(dose\_dicom\_dirs):
> print(f"Loading Dose from: {os.path.basename(dose\_dir)}”)
> dose\_image, ref\_file = read\_dicom\_series(dose\_dir, is\_dose=True)
> dose\_images.append(dose\_image)
> dose\_ref\_files.append(ref\_file)
> print(f" Dose size: {dose\_image.GetSize()}“)
> # Selects reference image
> reference\_ct, ref\_index = select\_reference\_image(ct\_images, ct\_names)
> print(”\n— Performing Registration & Resampling —“)
> # Process each dataset
> for i in range(len(ct\_images)):
> print(f”\n{‘-’\*50}“)
> print(f"Processing dataset {i+1}/{len(ct\_images)}: {ct\_names[i]}”)
> # Gets transform
> if i == ref\_index:
> transform = sitk.Euler3DTransform() # Identity transform for reference
> print(" Using identity transform (reference image)“)
> else:
> print(f” Registering to reference image…“)
> transform = robust\_rigid\_registration(reference\_ct, ct\_images[i])
> # Setups resampler
> resampler = sitk.ResampleImageFilter()
> resampler.SetReferenceImage(reference\_ct)
> resampler.SetInterpolator(sitk.sitkLinear)
> resampler.SetTransform(transform)
> # Resamples CT
> print(” Resampling CT image…“)
> resampler.SetDefaultPixelValue(-1024.0)
> resampled\_ct = resampler.Execute(sitk.Cast(ct\_images[i], sitk.sitkFloat32))
> # Resamples dose
> print(” Resampling Dose image…“)
> resampler.SetDefaultPixelValue(0.0)
> resampled\_dose = resampler.Execute(sitk.Cast(dose\_images[i], sitk.sitkFloat64))
> # Creates body mask
> print(” Creating body mask…“)
> body\_mask = create\_body\_mask(resampled\_ct)
> # Creates output directory for this dataset
> dataset\_output\_dir = os.path.join(output\_dir, ct\_names[i])
> os.makedirs(dataset\_output\_dir, exist\_ok=True)
> # Saves as NIfTI
> print(” Saving NIfTI files…“)
> sitk.WriteImage(resampled\_ct, os.path.join(dataset\_output\_dir, f’{ct\_names[i]}\_CT.nii.gz’))
> sitk.WriteImage(resampled\_dose, os.path.join(dataset\_output\_dir, f’{ct\_names[i]}\_dose.nii.gz’))
> sitk.WriteImage(body\_mask, os.path.join(dataset\_output\_dir, f’{ct\_names[i]}\_mask.nii.gz’))
> print(f” Saved: {ct\_names[i]}\_CT.nii.gz")
> print(f" Saved: {ct\_names[i]}\_dose.nii.gz")
> print(f" Saved: {ct\_names[i]}\_mask.nii.gz")
> # Process RT Structure
> rt\_struct\_dir = rt\_struct\_dirs[i]
> if os.path.exists(rt\_struct\_dir):
> print(f" Processing RT Structure from: {os.path.basename(rt\_struct\_dir)}“)
> rt\_files = [f for f in os.listdir(rt\_struct\_dir) if f.lower().endswith(‘.dcm’)]
> if rt\_files:
> rt\_file\_path = os.path.join(rt\_struct\_dir, rt\_files[0])
> success = transform\_rt\_structure(rt\_file\_path, transform, reference\_ct,
> dataset\_output\_dir, ct\_names[i])
> if success:
> print(f” Saved: {ct\_names[i]}\_RT\_Structure.dcm")
> else:
> print(f" Warning: No DICOM files found in {rt\_struct\_dir}“)
> else:
> print(f” Warning: RT Structure directory not found: {rt\_struct\_dir}“)
> print(f”\n{‘-’\*50}“)
> print(”=== WORKFLOW COMPLETED SUCCESSFULLY ===“)
> print(f"Output directory: {output\_dir}”)
> print(“All files saved in NIfTI format (.nii.gz)”)
> if **name** == “ **main** ”:
> main()
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [August 25, 2025, 5:11pm UTC](https://discourse.itk.org/t/unifying-ct-images-geometry-and-cropping-to-common-roi/7614/4 "2025-08-25T17:11:31Z")
Hello @blue_sky,
Based on the description it appears that the problem is with saving the non-image, RTStruct, DICOM data. SimpleITK supports operations and IO for images and you should be able to correctly save the results in DICOM format (see [this example](https://simpleitk.readthedocs.io/en/master/link_DicomSeriesReadModifyWrite_docs.html)).
Saving the RTStruct is not supported by SimpleITK, so possibly use pydicom (see [this discussion](https://github.com/pydicom/pydicom/issues/1122)), dedicated tools such as [RT-Utils](https://github.com/qurit/RT-Utils) or Slicer’s [Slicer RT extension](https://slicerrt.github.io/) which enables this form of [DICOM export](https://slicer.readthedocs.io/en/latest/user_guide/modules/segmentations.html#dicom-export).