Resampling two volumes with different number of slices.

DushiFdo · June 3, 2022, 8:29pm

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:

Should I register the two phases first and then resample? or vice versa?
Why some slices are full black in the resampled image?
I am planning to convert to isotrophic voxel spacing later (after resampling and registration). Is that correct?

zivy · June 4, 2022, 7:18pm

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:

Register original images.
Define the isotropic grid which will be associated with the venous image (origin, size, isotropic-spacing, direction cosine matrix).
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:

Fundamental concepts.
Registration.
Skim through the registration notebooks in the Jupyter notebook repository (Python notebooks the 6* series of notebooks).
Function for making images isotropic, available in this notebook (called make_isotropic).
If you have the time, go over the SimpleITK tutorial.

DushiFdo · June 5, 2022, 3:38am

Hi @zivy ,

Thank you so much for the detailed explanation. I registered two volumes following registration notebooks as follows:

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'))

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:

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?

zivy · June 5, 2022, 3:12pm

Hello @DushiFdo,

Let’s start with the transformation, it isn’t what you think it is (see this notebook). 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. 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. See the program’s documentation on registration.

DushiFdo · June 5, 2022, 6:08pm

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:

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:

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.

Why/how to avoid registered image ending up with gray pixels in some parts of the image?
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)?
Is there a way to prevent the resampled image to same size (512,512,156) prevent from ending up with full black slices?

zivy · June 6, 2022, 1:11pm

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, code:

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 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.

DushiFdo · June 6, 2022, 3:58pm

@zivy ,

Thanks again for the reply. Based on you recommendations, these are the steps I followed:

Register original images

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.

(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)

Next I transformed the fixed image so that it has isotropic voxel values:

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:

(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)

I then resampled the moving_resampled image (obtained after registration) to isotopic fixed image domain:

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. ?

zivy · June 6, 2022, 8:57pm

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:

isotropic_moving = sitk.Resample(moving_image, isotropic_fixed_image, final_transform)

DushiFdo · June 7, 2022, 12:36am

Hi @zivy,

The above code still does not work (corresponding slices in the two phases do not match after incorporating final_transform from registration).

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:

jizhang02 · January 11, 2023, 10:49am

Hello,
I have the same task as you, which is resampling, I would like to know if you have solved this problem. Thank you!

DushiFdo · January 11, 2023, 4:45pm

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)