# Resampling two volumes with different number of slices.
**URL:** https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051
**Category:** Beginner Questions
**Tags:** registration, dicom, simpleitk
**Created:** [June 3, 2022, 8:29pm UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051 "2022-06-03T20:29:43Z")
**Posts on this page:** 11
**Page:** 1
### Author:  [@DushiFdo](https://discourse.itk.org/u/DushiFdo)
#### Post date: [June 3, 2022, 8:29pm UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051/1 "2022-06-03T20:29:43Z")
I have two CT volumetric scans for the same patient (for arterial and venous) with different number of slices. e.g. arterial scan has ~146 dicom images whereas venous scan has ~156 dicom images.
I want to resample one phase (e.g. arterial) to another (e.eg venous) so that they have the same number of images and slices are aligned (for example, slice number 100 in arterial looks same as that of slice number 100 in venous).
I tried the following:
`art_image_resampled = sitk.Resample(art_image, ven_image)` but some images in arterial phase are full black.
My questions are as follows:
1. Should I register the two phases first and then resample? or vice versa?
2. Why some slices are full black in the resampled image?
3. I am planning to convert to isotrophic voxel spacing later (after resampling and registration). Is that correct?
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [June 4, 2022, 7:18pm UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051/2 "2022-06-04T19:18:05Z")
Hello @DushiFdo,
Welcome to SimpleITK!
The reason some slices in the arterial phase are black is that those are voxels from the venous phase that are falling outside the arterial phase CT when the mapping between the two is the identity transform.
Recommended workflow:
1. Register original images.
2. Define the isotropic grid which will be associated with the venous image (origin, size, isotropic-spacing, direction cosine matrix).
3. Use the grid settings to resample the venous CT with the identity transform. Use the grid settings to resample the arterial CT with the transform obtained from the registration.
Some recommended reading:
1. [Fundamental concepts](https://simpleitk.readthedocs.io/en/master/fundamentalConcepts.html).
2. [Registration](https://simpleitk.readthedocs.io/en/master/registrationOverview.html).
3. Skim through the registration notebooks in the [Jupyter notebook repository](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks) (Python notebooks the `6*` series of notebooks).
4. Function for making images isotropic, available in [this notebook](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/master/Python/05_Results_Visualization.ipynb) (called `make_isotropic`).
5. If you have the time, go over the [SimpleITK tutorial](https://simpleitk.org/TUTORIAL/).
---
### Author:  [@DushiFdo](https://discourse.itk.org/u/DushiFdo)
#### Post date: [June 5, 2022, 3:38am UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051/3 "2022-06-05T03:38:46Z")
Hi @zivy ,
Thank you so much for the detailed explanation. I registered two volumes following registration notebooks as follows:
```python
data_directory = r"venous"
series_IDs = sitk.ImageSeriesReader.GetGDCMSeriesIDs(data_directory)
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(data_directory, series_IDs[0])
data_directory2 = r"arterial"
series_IDs2 = sitk.ImageSeriesReader.GetGDCMSeriesIDs(data_directory2)
series_file_names2 = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(data_directory2, series_IDs2[0])
fixed_image = sitk.ReadImage(series_file_names, sitk.sitkFloat32)
moving_image = sitk.ReadImage(series_file_names2, sitk.sitkFloat32)
initial_transform = sitk.CenteredTransformInitializer(fixed_image,
moving_image,
sitk.ScaleVersor3DTransform(),
sitk.CenteredTransformInitializerFilter.GEOMETRY)
moving_resampled = sitk.Resample(moving_image, fixed_image, initial_transform, sitk.sitkLinear, 0.0, moving_image.GetPixelID())
registration_method = sitk.ImageRegistrationMethod()
# Similarity metric settings.
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.01)
registration_method.SetInterpolator(sitk.sitkLinear)
# Optimizer settings.
registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, convergenceMinimumValue=1e-6,
convergenceWindowSize=10)
registration_method.SetOptimizerScalesFromPhysicalShift()
# Don't optimize in-place, we would possibly like to run this cell multiple times.
registration_method.SetInitialTransform(initial_transform, inPlace=False)
# Connect all of the observers so that we can perform plotting during registration.
registration_method.AddCommand(sitk.sitkStartEvent, start_plot)
registration_method.AddCommand(sitk.sitkEndEvent, end_plot)
registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, update_multires_iterations)
registration_method.AddCommand(sitk.sitkIterationEvent, lambda: plot_values(registration_method))
final_transform = registration_method.Execute(sitk.Cast(fixed_image, sitk.sitkFloat32),
sitk.Cast(moving_image, sitk.sitkFloat32))
print(f'Final metric value: {registration_method.GetMetricValue()}')
print(f'Optimizer\'s stopping condition, {registration_method.GetOptimizerStopConditionDescription()}')
moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform, sitk.sitkLinear, 0.0, moving_image.GetPixelID())
interact(display_images, fixed_image_z=(0,fixed_image.GetSize()[2]-1), moving_image_z=(0,moving_image.GetSize()[2]-1),
fixed_npa = fixed(sitk.GetArrayViewFromImage(fixed_image)), moving_npa=fixed(sitk.GetArrayViewFromImage(moving_image)));
interact(display_images_with_alpha, image_z=(0,fixed_image.GetSize()[2] - 1), alpha=(0.0,1.0,0.05),
fixed = fixed(fixed_image), moving=fixed(moving_resampled));
sitk.WriteImage(moving_resampled, os.path.join(OUTPUT_DIR, 'arterial_resampled.mha'))
#print(moving_resampled)
sitk.WriteTransform(final_transform, os.path.join(OUTPUT_DIR, 'venous_2_arterial.tfm'))
```
1. Here, is my co-registered image now is given by ‘‘arterial\_resampled.mha’’? While the resampled image has the same dimensions as that of venous, it gives weird looking images as follows:
My code to resample dicom images to istropic voxel sizes as follows:
```auto
def resample(image, scan, new_spacing=[1,1,1]):
spacing = np.array([float(scan[0].SliceThickness),
float(scan[0].PixelSpacing[0]),
float(scan[0].PixelSpacing[0])])
resize_factor = spacing / new_spacing
new_real_shape = image.shape * resize_factor
new_shape = np.round(new_real_shape)
real_resize_factor = new_shape / image.shape
new_spacing = spacing / real_resize_factor
image = scipy.ndimage.interpolation.zoom(image, real_resize_factor)
return image
```
Could you please explain better/provide an example for points 2 and 3 in your answer above?
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [June 5, 2022, 3:12pm UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051/4 "2022-06-05T15:12:11Z")
Hello @DushiFdo,
Let’s start with the transformation, it isn’t what you think it is (see [this notebook](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/master/Python/22_Transforms.ipynb)). I suspect you want to use `ComposeScaleSkewVersor3DTransform`. Similarity metric is fine, but possibly you can get by using correlation (`SetMetricAsCorrelation`) which is computationally less expensive.
To evaluate registration results I prefer to use a linked cursor approach, save the resampled image to disk and open both images in [ITK-SNAP](http://www.itksnap.org/pmwiki/pmwiki.php). When you click on one point in one instance of the program it will jump to the same spatial location in the other. You can then judge if the points in both match up.
Your code for resampling is incorrect. The SliceThickness DICOM tag is not the slice spacing which is what should be used. I highly recommend using the function I pointed to which takes care of other subtleties such as a non-identity direction cosine matrix. Original CT images may not be axis aligned, a scipy assumption.
Given your struggles, possibly consider using [Slicer](https://www.slicer.org/). See the program’s documentation on [registration](https://slicer.readthedocs.io/en/latest/user_guide/registration.html).
---
### Author:  [@DushiFdo](https://discourse.itk.org/u/DushiFdo)
#### Post date: [June 5, 2022, 6:08pm UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051/5 "2022-06-05T18:08:29Z")
Hi @zivy,
I tried registration with `ComposeScaleSkewVersor3DTransform` (I don’t totally understand how to choose the best transformation) but I still get gray pixels in the registered image as follows:
As per your advice, I resample to isotrophic images as follows:
```python
def resample_image(image, out_spacing=(1.0, 1.0, 1.0), out_size=None, is_label=False, pad_value=0):
"""Resamples an image to given element spacing and output size."""
original_spacing = np.array(image.GetSpacing())
original_size = np.array(image.GetSize())
if out_size is None:
out_size = np.round(np.array(original_size * original_spacing / np.array(out_spacing))).astype(int)
else:
out_size = np.array(out_size)
original_direction = np.array(image.GetDirection()).reshape(len(original_spacing),-1)
original_center = (np.array(original_size, dtype=float) - 1.0) / 2.0 * original_spacing
out_center = (np.array(out_size, dtype=float) - 1.0) / 2.0 * np.array(out_spacing)
original_center = np.matmul(original_direction, original_center)
out_center = np.matmul(original_direction, out_center)
out_origin = np.array(image.GetOrigin()) + (original_center - out_center)
resample = sitk.ResampleImageFilter()
resample.SetOutputSpacing(out_spacing)
resample.SetSize(out_size.tolist())
resample.SetOutputDirection(image.GetDirection())
resample.SetOutputOrigin(out_origin.tolist())
resample.SetTransform(sitk.Transform())
resample.SetDefaultPixelValue(pad_value)
if is_label:
resample.SetInterpolator(sitk.sitkNearestNeighbor)
else:
resample.SetInterpolator(sitk.sitkLinear)
return resample.Execute(sitk.Cast(image, sitk.sitkFloat32))
```
Also, I tried re-sampling the arterial image(512,512,146) to venous image size(512,512,156) using the following:
```auto
def scale_resample(image, interpolator, new_size, new_spacing):
original_physical_size = [(sz-1)*spc for sz,spc in zip(image.GetSize(), image.GetSpacing())]
new_physical_size = [(sz-1)*spc for sz,spc in zip(new_size, new_spacing)]
scale_tx = sitk.ScaleTransform(image.GetDimension(), [ops/nps for nps, ops in zip(new_physical_size, original_physical_size)])
rigid_tx = sitk.Euler3DTransform()
rigid_tx.SetMatrix(image.GetDirection())
rigid_tx.SetTranslation(image.GetOrigin())
#accomodate for non identity direction cosine matrix, ensure that scaling is
#along the image axes and not the canonical ones.
tx = sitk.CompositeTransform([rigid_tx, scale_tx, rigid_tx.GetInverse()])
#output image has the same origin and direction as the input, but is a
#scaled version of the content
return sitk.Resample(image,
size = new_size,
transform = tx,
interpolator = interpolator,
outputOrigin = image.GetOrigin(),
outputDirection = image.GetDirection(),
outputSpacing = new_spacing)
```
But I still get black images for all the slices up-to slice 60.
1. Why/how to avoid registered image ending up with gray pixels in some parts of the image?
2. Some slices are full black after the co-registration. Is it because the depth in the fixed image(156) is bigger than the moving image(146)?
3. Is there a way to prevent the resampled image to same size (512,512,156) prevent from ending up with full black slices?
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [June 6, 2022, 1:11pm UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051/6 "2022-06-06T13:11:54Z")
Hello @DushiFdo,
You are thinking of registration+resampling+display as an atomic/single step operation. This is not the case.
**Registration** estimates a transformation mapping points between the fixed image to the moving image. Please see this [jupyter notebook](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/master/Python/63_Registration_Initialization.ipynb), code:
```auto
gui.RegistrationPointDataAquisition(
fixed_image,
moving_image,
figure_size=(8, 4),
known_transformation=final_transform,
fixed_window_level=ct_window_level,
moving_window_level=mr_window_level,
);
```
This code uses the computed transformation to map points between the two images without performing resampling. This way you can confirm that the transformation estimated via registration makes sense. If it does not then no need to do resampling as the registration has failed. If the mapping makes sense then you can resample.
**Resampling** , using the transformation from the previous step we resample the moving image onto the fixed image grid. Some points from the fixed image may map outside the moving image. The intensity value for these is set using a default value in the `Resample` function, often zero. When working with CT a more reasonable value is -1000, the Hounsfield value for air.
**Display** , we shouldn’t refer to the scalar intensity values based on their visualized color (black/gray) because that is what we see after the high dynamic range of the [Hounsfield scale](https://en.wikipedia.org/wiki/Hounsfield_scale) is mapped to a low dynamic range of [0,255] for display.
So, (1,2) gray pixels and “full black slices” are likely because the estimated transformation maps those points outside the moving image. Likely the intensity value for both is 0, but in slices that contain information it visually appears as gray and in those that are completely empty it appears as black, even though it is the same intensity value (remember the display step). (3) If the registration transformation maps all slice points completely outside the moving image, then no. Likely something is not working well with the registration step.
---
### Author:  [@DushiFdo](https://discourse.itk.org/u/DushiFdo)
#### Post date: [June 6, 2022, 3:58pm UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051/7 "2022-06-06T15:58:49Z")
@zivy ,
Thanks again for the reply. Based on you recommendations, these are the steps I followed:
1. Register original images
```python
data_directory = r"venous"
series_IDs = sitk.ImageSeriesReader.GetGDCMSeriesIDs(data_directory)
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(data_directory, series_IDs[0])
data_directory2 = r"arterial"
series_IDs2 = sitk.ImageSeriesReader.GetGDCMSeriesIDs(data_directory2)
series_file_names2 = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(data_directory2, series_IDs2[0])
fixed_image = sitk.ReadImage(series_file_names, sitk.sitkFloat32)
moving_image = sitk.ReadImage(series_file_names2, sitk.sitkFloat32)
initial_transform = sitk.CenteredTransformInitializer(fixed_image,
moving_image,
sitk.ComposeScaleSkewVersor3DTransform(),
sitk.CenteredTransformInitializerFilter.GEOMETRY)
moving_resampled = sitk.Resample(moving_image, fixed_image, initial_transform, sitk.sitkLinear, -1000, moving_image.GetPixelID())
registration_method = sitk.ImageRegistrationMethod()
# Similarity metric settings.
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.01)
registration_method.SetInterpolator(sitk.sitkLinear)
# Optimizer settings.
registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, convergenceMinimumValue=1e-6,
convergenceWindowSize=10)
registration_method.SetOptimizerScalesFromPhysicalShift()
# Don't optimize in-place, we would possibly like to run this cell multiple times.
registration_method.SetInitialTransform(initial_transform, inPlace=False)
# Connect all of the observers so that we can perform plotting during registration.
registration_method.AddCommand(sitk.sitkStartEvent, start_plot)
registration_method.AddCommand(sitk.sitkEndEvent, end_plot)
registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, update_multires_iterations)
registration_method.AddCommand(sitk.sitkIterationEvent, lambda: plot_values(registration_method))
final_transform = registration_method.Execute(sitk.Cast(fixed_image, sitk.sitkFloat32),
sitk.Cast(moving_image, sitk.sitkFloat32))
print(f'Final metric value: {registration_method.GetMetricValue()}')
print(f'Optimizer\'s stopping condition, {registration_method.GetOptimizerStopConditionDescription()}')
moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform, sitk.sitkLinear, 0.0, moving_image.GetPixelID())
interact(display_images, fixed_image_z=(0,fixed_image.GetSize()[2]-1), moving_image_z=(0,moving_image.GetSize()[2]-1),
fixed_npa = fixed(sitk.GetArrayViewFromImage(fixed_image)), moving_npa=fixed(sitk.GetArrayViewFromImage(moving_image)));
interact(display_images_with_alpha, image_z=(0,fixed_image.GetSize()[2] - 1), alpha=(0.0,1.0,0.05),
fixed = fixed(fixed_image), moving=fixed(moving_resampled));
```
After registration, the two images are aligned and have same domain: size, spacing, orgin, direction.
```auto
(512, 512, 156)
(0.869140625, 0.869140625, 3.0)
(-231.0654296875, -391.0654296875, 37.5)
(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
```
1. Next I transformed the fixed image so that it has isotropic voxel values:
```python
def resample_image(image, out_spacing=(1.0, 1.0, 1.0), out_size=None, is_label=False, pad_value=0):
"""Resamples an image to given element spacing and output size."""
original_spacing = np.array(image.GetSpacing())
original_size = np.array(image.GetSize())
if out_size is None:
out_size = np.round(np.array(original_size * original_spacing / np.array(out_spacing))).astype(int)
else:
out_size = np.array(out_size)
original_direction = np.array(image.GetDirection()).reshape(len(original_spacing),-1)
original_center = (np.array(original_size, dtype=float) - 1.0) / 2.0 * original_spacing
out_center = (np.array(out_size, dtype=float) - 1.0) / 2.0 * np.array(out_spacing)
original_center = np.matmul(original_direction, original_center)
out_center = np.matmul(original_direction, out_center)
out_origin = np.array(image.GetOrigin()) + (original_center - out_center)
resample = sitk.ResampleImageFilter()
resample.SetOutputSpacing(out_spacing)
resample.SetSize(out_size.tolist())
resample.SetOutputDirection(image.GetDirection())
resample.SetOutputOrigin(out_origin.tolist())
resample.SetTransform(sitk.Transform())
resample.SetDefaultPixelValue(pad_value)
if is_label:
resample.SetInterpolator(sitk.sitkNearestNeighbor)
else:
#resample.SetInterpolator(sitk.sitkBSpline)
resample.SetInterpolator(sitk.sitkLinear)
return resample.Execute(sitk.Cast(image, sitk.sitkFloat32))
```
Now the image domain is as follows:
```auto
(445, 445, 468)
(1.0, 1.0, 1.0)
(-231.0, -391.0, 36.5)
(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
```
1. I then resampled the moving\_resampled image (obtained after registration) to isotopic fixed image domain:
```python
resnew = sitk.Resample(moving_resampled, isotropic_image)
```
While this image has same domain as the isotropic image, slices do not overlap. e.g: slice 100 in isotropic fixed image is different than that of slice 100 in resampled moving image onto istropic image.
However, if I do `sitk.Resample(moving_image, isotropic_image)` they do overlap.
What am I doing wrong here? Why does slices do not overlap after registration even if I resample to `isotropic_image`. ?
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [June 6, 2022, 8:57pm UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051/8 "2022-06-06T20:57:53Z")
Hello @DushiFdo,
Hard to follow as the code is pretty long. In any case, you do not need to compute the `moving_resampled` image immediately after registration. The output from registration is just `final_transform`.
The next step is to resample the fixed image to have isotropic spacing, and you get `isotropic_fixed_image` out of that and finally resample the original moving image:
```auto
isotropic_moving = sitk.Resample(moving_image, isotropic_fixed_image, final_transform)
```
---
### Author:  [@DushiFdo](https://discourse.itk.org/u/DushiFdo)
#### Post date: [June 7, 2022, 12:36am UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051/9 "2022-06-07T00:36:05Z")
Hi @zivy,
The above code still does not work (corresponding slices in the two phases do not match after incorporating `final_transform` from registration).
```python
isotropic_image = resample_image(fixed_image)
isotropic_moving = sitk.Resample(moving_image, isotropic_image, final_transform)
```
For e.g., slice 400 from each image looks as follows:
---
### Author:  [@jizhang02](https://discourse.itk.org/u/jizhang02)
#### Post date: [January 11, 2023, 10:49am UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051/10 "2023-01-11T10:49:24Z")
Hello,
I have the same task as you, which is resampling, I would like to know if you have solved this problem. Thank you!
---
### Author:  [@DushiFdo](https://discourse.itk.org/u/DushiFdo)
#### Post date: [January 11, 2023, 4:45pm UTC](https://discourse.itk.org/t/resampling-two-volumes-with-different-number-of-slices/5051/11 "2023-01-11T16:45:29Z")
Hi, I didn’t have to register the two phases because the two phases seemed to be in the same coordinate system. To align the phases, I did the following, which worked for me:
`img1_resampled = sitk.Resample(img1, img2, defaultPixelValue=-1000)`