Resampling MR image using sitk.ResampleImageFilter() and a composite transform results in grey artefacts in output image

Question

I performed with FSL and ANTs registrations between the slices of a interleaved aquired t2 weighted MR image to get rid of geometric distortions (rotation, translation..) and i wanted to compute a composite transformation using sitk composite transform for this approach. I inspected the resulting Composite transform and it should work fine, but i have in both methods (FSL and ANTs) the issue that the resampled image has grey artefacts in the last slices or in parts of the last slice.Regardless if i only add one transformation to the composite transform or 40. I add the identity matrix, then the matrix from slice one and two and then all transformation matrices from each slice sequentially.

The SetUseNearestNeighborExtrapolator() from SimpleITK does adress this problem, but as a result it only maps the pixels from my reference image to the registered image where the moving image pixels fall out of the grid.

I dont know what to try anymore because i tried differend composite transformation computation approaches but i always have some of these artefacts. The impact of the artefact is different for different images.

def create_transformations(even_matrices, odd_matrices, matrix_12):
            transformations = []
            # Add initial rotation and translation

            # Create an identity transform
            identity_matrix = np.eye(4)
            identity_affine_transform = sitk.AffineTransform(3)
            identity_affine_transform.SetMatrix(identity_matrix[:3, :3].flatten())
            identity_affine_transform.SetTranslation(identity_matrix[:3, 3])
            sitk.WriteTransform(identity_affine_transform, os.path.join(output_directory, 'identity_transform.txt'))
            transformations.append(identity_affine_transform)

            
            initial_rotation_matrix = np.array(matrix_12.GetMatrix()).reshape((3, 3)).T
            initial_translation_vector = np.array(matrix_12.GetTranslation())
            initial_affine_transform = sitk.AffineTransform(3)
            initial_affine_transform.SetMatrix(initial_rotation_matrix.flatten())
            initial_affine_transform.SetTranslation(initial_translation_vector)
            transformations.append(initial_affine_transform)


            max_len = max(len(even_matrices), len(odd_matrices))

            for i in range(max_len):
                if i < len(odd_matrices):
                    odd_transform = odd_matrices[i]
                    
                    # Convert the affine transform to a homogeneous matrix
                    homogeneous_matrix = get_affine_homogeneous_matrix(3, odd_transform)
                    
                    # Extract the matrix and translation from the homogeneous matrix
                    odd_matrix, translation = get_affine_matrix_and_translation(3, homogeneous_matrix)
                    
                    # Create and append the transform again to ensure it's correctly set
                    odd_affine_transform = sitk.AffineTransform(3)
                    odd_affine_transform.SetMatrix(odd_matrix)
                    odd_affine_transform.SetTranslation(translation)
                    sitk.WriteTransform(odd_affine_transform, os.path.join(output_directory,  f"{i}_odd_transform.txt"))
                    transformations.append(odd_affine_transform)

                if i < len(even_matrices):
                    even_transform = even_matrices[i]
                    
                    # Convert the affine transform to a homogeneous matrix
                    homogeneous_matrix = get_affine_homogeneous_matrix(3, even_transform)
                    
                    # Extract the matrix and translation from the homogeneous matrix
                    even_matrix, translation = get_affine_matrix_and_translation(3, homogeneous_matrix)
                    
                    # Create and append the transform again to ensure it's correctly set
                    even_affine_transform = sitk.AffineTransform(3)
                    even_affine_transform.SetMatrix(even_matrix)
                    even_affine_transform.SetTranslation(translation)
                    sitk.WriteTransform(even_affine_transform, os.path.join(output_directory,  f"{i}_even_transform.txt"))
                    transformations.append(even_affine_transform)
                     # Process all odd matrices first
            
           

            return transformations

def get_affine_homogeneous_matrix(dim, transformation):
    """
    Returns a homogeneous matrix for an affine transformation.
    :param dim: The dimension of the transformation.
    :param transformation: The sitk transformation.
    :return: A homogeneous (dim+1)x(dim+1) matrix as an np.array.
    """
    matrix = np.eye(dim + 1)
    matrix[:dim, :dim] = np.array(transformation.GetMatrix()).reshape([dim, dim]).T
    matrix[dim, :dim] = np.array(transformation.GetTranslation())
    return matrix


def get_affine_matrix_and_translation(dim, homogeneous_matrix):
    """
    Returns an affine transformation parameters for a homogeneous matrix.
    :param dim: The dimension of the transformation.
    :param homogeneous_matrix: The homogeneous (dim+1)x(dim+1) matrix as an np.array.
    :return: A tuple of the homogeneous matrix as a list, and the translation parameters as a list.
    """
    matrix = homogeneous_matrix[:dim, :dim].T.reshape(-1).tolist()
    translation = homogeneous_matrix[dim, :dim].tolist()
    return matrix, translation


def create_composite(dim, transformations, merge_affine=False):
    """
    Creates a composite sitk transform based on a list of sitk transforms.
    :param dim: The dimension of the transformation.
    :param transformations: A list of sitk transforms.
    :param merge_affine: If true, merge affine transformations before calculating the composite transformation.
    :return: The composite sitk transform.
    """
    if merge_affine:
        merged_transformations = []
        combined_matrix = None
        for transformation in transformations:
            if isinstance(transformation, sitk.AffineTransform):
                if combined_matrix is None:
                    combined_matrix = np.eye(dim + 1)
                current_matrix = get_affine_homogeneous_matrix(dim, transformation)
                combined_matrix = current_matrix @ combined_matrix
            else:
                if combined_matrix is not None:
                    matrix, translation = get_affine_matrix_and_translation(dim, combined_matrix)
                    combined_affine_transform = sitk.AffineTransform(dim)
                    combined_affine_transform.SetMatrix(matrix)
                    combined_affine_transform.SetTranslation(translation)
                    merged_transformations.append(combined_affine_transform)
                merged_transformations.append(transformation)
                combined_matrix = None
        if combined_matrix is not None:
            matrix, translation = get_affine_matrix_and_translation(dim, combined_matrix)
            combined_affine_transform = sitk.AffineTransform(dim)
            combined_affine_transform.SetMatrix(matrix)
            combined_affine_transform.SetTranslation(translation)
            merged_transformations.append(combined_affine_transform)
        transformations = merged_transformations

    if sitk.Version_MajorVersion() == 1:
        compos = sitk.Transform(dim, sitk.sitkIdentity)
        for transformation in transformations:
            compos.AddTransform(transformation)
    else:
        compos = sitk.CompositeTransform(transformations)
    return compos

def save_composite_transform(transformations, filename, output_directory):
    composite_transform = create_composite(3, transformations, merge_affine=True)
    sitk.WriteTransform(composite_transform, os.path.join(output_directory, filename))
    return composite_transform

This results in a composite transform which looks fine when i print it out but when i perform resampling the image do consist always grey artefacts.

Answer 1

It is hard to identify what is not working, given the lengthy code, and the specific issue "grey artefacts". Educated guess is that these are due to samples mapped outside the image and this is the default value returned. Why they are mapped outside the image is unclear if the transform is correct.

I would suggest asking the question on the dedicated ITK discourse forum after reducing the code to a minimal working example and including a sample image illustrating the problem.

An additional comment, it appears the code is also supporting SimpleITK versions from the first major release 1.x. These are ancient (2.0.0 was released four years ago) and unsupported versions, highly recommend not using them.