3D Affine registration: different datasets

[DushiFdo (Dushi)]

Hello,

I am working on 3D CT images. I need to perform registration on two different datasets. I randomly selected a scan from dataset 1 as the fixed_image and perform a 3d affine registration to register all the scans from dataset 2 to it. My code is as follows:

def register(fixed_image, moving_image):
    initial_transform = sitk.CenteredTransformInitializer(fixed_image,
                                                      moving_image,
                                                      sitk.AffineTransform(fixed_image.GetDimension()),
                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)

    moving_resampled = sitk.Resample(moving_image, fixed_image, initial_transform, sitk.sitkLinear, -1000, moving_image.GetPixelID())
    
    registration_method = sitk.ImageRegistrationMethod()
 
    registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
    registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
    registration_method.SetMetricSamplingPercentage(0.01)

    registration_method.SetInterpolator(sitk.sitkLinear)

    registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=500, convergenceMinimumValue=1e-6, 
                                                       convergenceWindowSize=10)
    registration_method.SetOptimizerScalesFromPhysicalShift()
                                     
    registration_method.SetInitialTransform(initial_transform, inPlace=False)

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

    moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform, sitk.sitkLinear,-1000,moving_image.GetPixelID())                         
    moving_resampled_npa = sitk.GetArrayFromImage(moving_resampled)
   
    return moving_resampled_npa   

However, after registration there is a gray color border around some images. How can I avoid it?

In the above image, first row is after registration. Second row is before registration.

[zivy (Ziv Yaniv)]

Hello @DushiFdo,

This is the default value you selected for intensities that fall outside the resampled volume (-1000, the HU value for air which is specified in the call to Resample). You can see that the intensity matches the air around the patient. If you want it to match the border surrounding the original image you need to use that value instead.

[DushiFdo (Dushi)]

Hi @zivy ,

This is how a slice from fixed image looks like:

I changed the -1000 value to 0.0 as follows:

moving_resampled = sitk.Resample(moving_image, fixed_image, initial_transform, sitk.sitkLinear, 0.0, moving_image.GetPixelID())

and,

moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform, sitk.sitkLinear,0.0,moving_image.GetPixelID())

Still the resulting image has gray borders (first row):

[zivy (Ziv Yaniv)]

Hello @DushiFdo,

Not sure if 0.0 is the right value to get the effect you want. In the image below the region marked in green has a value of -1000, it’s air, the region marked in red has another value which is likely not 0, but is the value you want. You need to use that value.

[DushiFdo (Dushi)]

Thanks @zivy . Is there a way to extract that specific value?

[zivy (Ziv Yaniv)]

Open the file in ITK-SNAP, place cursor at desired location and on the left side you will see the “Intensity under cursor”.

[DushiFdo (Dushi)]

@Zivy,

Intensity around the image (black area) has different values at different points like -997, -103, -1024 and so on. What can I do in this case?

[lassoan (Andras Lasso)]

In 3D Slicer, we use the median intensity value of the 8 corner voxels of the volume as background value and it works very well.

Slicer/Slicer/blob/a5f75351073ef62fd6198d9480d86c0009d70f9b/Libs/MRML/Core/vtkMRMLVolumeNode.cxx#L1167-L1199

[DushiFdo (Dushi)]

Thank you. How can I extract corner voxel values in python?

[lassoan (Andras Lasso)]

The Python API is the same as C++, there are just a few minor syntax differences (e.g., replace -> by .). So, you can use GetScalarComponentAsDouble method in Python as in the C++ example above.

[DushiFdo (Dushi)]

Thank you so much for your reply. This is the function I converted to python:

def GetImageBackgroundScalarComponentAsDouble(component):
    imageData = component.GetImageData()
    if imageData is None or component >= imageData.GetNumberOfScalarComponents():
        return 0.0

    extent = [0, -1, 0, -1, 0, -1]
    imageData.GetExtent(extent)

    if extent[0] > extent[1] or extent[2] > extent[3] or extent[4] > extent[5]:
        return 0.0

    scalarValues = []
    scalarValues.append(imageData.GetScalarComponentAsDouble(extent[0], extent[2], extent[4], component))
    scalarValues.append(imageData.GetScalarComponentAsDouble(extent[0], extent[2], extent[5], component))
    scalarValues.append(imageData.GetScalarComponentAsDouble(extent[0], extent[3], extent[4], component))
    scalarValues.append(imageData.GetScalarComponentAsDouble(extent[0], extent[3], extent[5], component))
    scalarValues.append(imageData.GetScalarComponentAsDouble(extent[1], extent[2], extent[4], component))
    scalarValues.append(imageData.GetScalarComponentAsDouble(extent[1], extent[2], extent[5], component))
    scalarValues.append(imageData.GetScalarComponentAsDouble(extent[1], extent[3], extent[4], component))
    scalarValues.append(imageData.GetScalarComponentAsDouble(extent[1], extent[3], extent[5], component))

    # Get the median value (nth_element performs partial sorting until nth largest element is found)
    MEDIANELEMENTINDEX = 3
    std::nth_element(scalarValues.begin(), scalarValues.begin() + MEDIANELEMENTINDEX, scalarValues.end())
    medianValue = scalarValues[MEDIANELEMENTINDEX]

    return medianValue

But I get a SyntaxError: invalid syntax :
std::nth_element(scalarValues.begin(), scalarValues.begin() + MEDIANELEMENTINDEX, scalarValues.end())

I am not sure what I am doing wrong here.

Also, should I call the function as:

GetImageBackgroundScalarComponentAsDouble(image_sitk)

[lassoan (Andras Lasso)]

std:: refers to the C++ standard template library, which is not available in Python (it is not related to VTK). Instead, you can use np.median function to get the median of the values in the scalarValues list.

[DushiFdo (Dushi)]

I changed it to np.median(scalarValues.begin(), scalarValues.begin() + MEDIANELEMENTINDEX, scalarValues.end()) and I called the function as: GetImageBackgroundScalarComponentAsDouble(image_sitk). My image is an sitk image and now I get, AttributeError: 'Image' object has no attribute 'GetImageData'

[lassoan (Andras Lasso)]

It seems that your image is not a VTK image but an ITK image - sorry, I could have thought, since it is the ITK forum after all.

So, you need to use ITK functions to get the voxel value. I’m not sure about the exact syntax, it is something like image.GetPixel((2,3,4)) - but it may depend on if you use SimpleITK or ITK-Python; and you may also choose get the voxels as a numpy array and use numpy indexing to get corner values.

[DushiFdo (Dushi)]

Thanks @lassoan .
@zivy, could you please tell me how to get the 8 corner voxel values in sitk so that I can get the median and set that as the default value for background as @lassoan suggested ? I tried to find it in the documentation, but couldn’t find any information.

[DushiFdo (Dushi)]

Hi @zivy, I tried normalizing the fixed image (so that the background voxel is 0.0) and set 0 as the default value, but I still get a gray border around the registered image. What could be the reason? Both fixed and moving image has the same dimensions (128x128x128).

[zivy (Ziv Yaniv)]

Hello @DushiFdo,

Image size/dimension has nothing to do with what you are trying to accomplish, not sure why you are providing that information. With respect to computing the median intensity value from the corners of an image (irrespective of dimensionality, so 4 in 2D, 8 in 3D), code below:

import SimpleITK as sitk
import numpy as np
import itertools

image = sitk.ReadImage(file_name)

corners_median_intensity = np.median([image[index] for index in itertools.product(*zip([0]*image.GetDimension(),[sz-1 for sz in image.GetSize()]))])

print(corners_median_intensity)

[DushiFdo (Dushi)]

Thank You. The median is -1024 (which I tried before as well). Still the image has gray borders (fixed and moving images are overlaid here):

Changes made (default value is set to -1024):

moving_resampled = sitk.Resample(moving_image, fixed_image, initial_transform, sitk.sitkLinear, -1024, moving_image.GetPixelID())

and,

moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform, sitk.sitkLinear,-1024,moving_image.GetPixelID())

I also tried changing the transformation from sitk.AffineTransform(fixed_image.GetDimension()) to
sitk.Similarity3DTransform() which resulted in narrower gray borders.

[mauigna06 (Mauro I. Dominguez)]

I think you need to sample the moving image corners median HU value before transforming, and use that value for the resample.

Hope it helps

[DushiFdo (Dushi)]

I figured that for images that looks like a gray circle around it, the background is -3024 even if the median of the 8 corner voxels is -1024.

For the images that looks as follows, it is -1024.

Is there a way to distinguish between the two without visualizing?