Understanding CompositeTransform discrepancies

Beginner Questions
,
1

I’m getting different output from the two code blocks below. The only difference is that in the second I wrap the rotation with a CompositeTransform. I would expect that this transform would just apply the rotation, but there is something weird going on that I don’t understand. Am I missing something?

Screenshot 2021-10-05 at 15.51.56
Screenshot 2021-10-05 at 15.51.561920×1099 143 KB

Screenshot 2021-10-05 at 15.52.06
Screenshot 2021-10-05 at 15.52.061920×1099 144 KB 

import SimpleITK as sitk
import matplotlib.pyplot as plt
import numpy as np

def myshow(img, title=None, margin=0.05, dpi=80):
    nda = sitk.GetArrayViewFromImage(img)[:,:,62]
    spacing = img.GetSpacing()
        
    ysize = nda.shape[0]
    xsize = nda.shape[1]
      
    figsize = ((1 + margin) * ysize / dpi, (1 + margin) * xsize / dpi)

    fig = plt.figure(title, figsize=figsize, dpi=dpi)
    ax = fig.add_axes([margin, margin, 1 - 2*margin, 1 - 2*margin])
    
    extent = (0, xsize*spacing[1], 0, ysize*spacing[0])

    t = ax.imshow(nda,
            extent=extent,
            interpolation='none',
            cmap='gray',
            origin='lower',
            vmin=0.0, vmax=1.0)
    
    if(title):
        plt.title(title)

    plt.show()

if __name__ == '__main__':
    scale = 1.0
    spacing = 0.2
    image_size = np.array([128, 250, 250])
    patch_size = np.array([64, 125, 125])
    sampling = sitk.sitkBSplineResamplerOrder3

    grid_array = ((np.indices(image_size) // 4).sum(axis=0) % 2).astype(float)
    grid = sitk.GetImageFromArray(grid_array)
    grid.SetSpacing(np.ones(3) * spacing)
    grid.SetOrigin((30.9, 10.2, 28.4))


    patch_border = 30
    patch_size_augmented = patch_size + 2*patch_border
    patch_origin = ((0.5 * image_size - 0.5 * patch_size_augmented) * np.ones(3)).astype(int)
    patch = grid_array[
        patch_origin[0]:patch_origin[0]+patch_size_augmented[0],
        patch_origin[1]:patch_origin[1]+patch_size_augmented[1],
        patch_origin[2]:patch_origin[2]+patch_size_augmented[2]
    ]
    patch_image = sitk.GetImageFromArray(patch)

    rotation = sitk.VersorTransform(np.array([1.0, 0.0, 0.0]), 0.15, 0.5 * np.array(patch.shape)[::-1])

    resampled_patch_image = sitk.Resample(patch_image, rotation, sampling, 255)
    resampled2 = sitk.GetArrayViewFromImage(resampled_patch_image)
    resampled2 = resampled2[
        patch_border:patch_border+patch_size[0],
        patch_border:patch_border+patch_size[1],
        patch_border:patch_border+patch_size[2]
    ]
    resampled2 = sitk.GetImageFromArray(resampled2)

    myshow(resampled2, 'Resampled')

    transform = []
    transform.append(rotation)
    transform = sitk.CompositeTransform(transform)
    resampled_patch_image = sitk.Resample(patch_image, transform, sampling, 255)
    resampled2 = sitk.GetArrayViewFromImage(resampled_patch_image)
    resampled2 = resampled2[
        patch_border:patch_border+patch_size[0],
        patch_border:patch_border+patch_size[1],
        patch_border:patch_border+patch_size[2]
    ]
    resampled2 = sitk.GetImageFromArray(resampled2)

    myshow(resampled2, 'Resampled Composite')
2

Hello @Theoretialperson,

You have uncovered a bug. The rotation is not cloned correctly inside the composite transform, the offset portion of the transformation is not updated.

Temporary solution: Use VersorRigid3DTransform with zero translation.

Minimal working example showing bug and temporary workaround:

    import SimpleITK as sitk
    
    axis = [1.0,0.0,0.0]
    angle = 0.15
    fixed_center = [92.5, 92.5, 62]

    # rotation is not copied correctly inside the composite transform
    # offset parameter is not updated
    rotation = sitk.VersorTransform(axis, angle, fixed_center)
    print(rotation)
    print('*'*30)
    print(sitk.CompositeTransform([rotation]))

    print('*'*30)

    # workaround: use a versor based rigid transformation with zero translation
    rigid = sitk.VersorRigid3DTransform(axis, angle, [0,0,0], fixed_center)
    print(rigid)
    print(sitk.CompositeTransform([rigid]))

Edit: bug reported here.

3 Likes