Transforming coordinates with saved affine matrix does not match transformed volume

[RuihanZhang2015]

I wish to transform coordinates (z,y,x) in one volume with two consecutive affine matrice. I also have a volume transformed with the same two affine matrices. However, I realized the transformed points do not appear in the position I expected. I wonder what is wrong with my code.

Here is the code for computing the affine matrix.

original_spacing = [0.4, 0.1625, 0.1625]
factors = (2, 4, 4)
new_spacing = original_spacing * np.array(factors)
fix_downsampled = downsample(fix_vol, factors)
move_downsampled = downsample(mov_vol, factors)      

ransac_kwargs = {'blob_sizes': [6,10]}

affine_kwargs = { 'metric' : 'MMI',
                        'optimizer':'LBFGSB',
                        'alignment_spacing': 1,
                        'shrink_factors': ( 4, 2, 1),
                        'smooth_sigmas': ( 0., 0., 0.),
                    }

steps = [('ransac', ransac_kwargs)]

affine = alignment_pipeline(fix_downsampled, move_downsampled, new_spacing, new_spacing, steps)
affine2 = alignment_pipeline(fix_vol, mov_vol, original_spacing, original_spacing, steps=[('affine', affine_kwargs)],static_transform_list=[affine])

os.makedirs(os.path.dirname(affine_filename),exist_ok=True)
with open(affine_filename,'wb') as f:
    pickle.dump([affine,affine2],f)
affine_list = [affine,affine2]

Here is the code for tranforming the volume.

affine_filename = f'{output_path}/affine_two_steps.pkl'
with open(affine_filename,'rb') as f:
        affine_list = pickle.load(f)

mov_h5_new = mov_h5.replace('_raw.h5','.h5')
print(mov_h5_new)
if os.path.exists(mov_h5_new):
    os.remove(mov_h5_new)

with h5py.File(mov_h5_new, "w") as f_new:
    for channel in ['405','488','561','594','640']:
        print(channel,flush = True)
        with h5py.File(mov_h5, "r") as f:
            print(f.keys())
            mov_vol = f[channel][:]
            aligned_vol = apply_transform(
                fix_vol, mov_vol,
                original_spacing, original_spacing,
                transform_list=affine_list,
            )
            if channel in f_new.keys():
                del f_new[channel]
            f_new.create_dataset(channel, aligned_vol.shape, dtype = aligned_vol.dtype, data=aligned_vol)
print(f'{mov_h5_new} finish transforming other channels',flush = True)  

Here is the code for transforming the coordinates:

spacing = [0.4, 0.1625, 0.1625]
original_spacing = np.array(spacing)

transformed_coords = apply_transform_to_coordinates(
    coordinates = original_coords,
    transform_list = affine_list,
    transform_spacing = original_spacing,
    transform_origin = None,
)

Here is the results. The expected points should appear in the following location:

transformed [[  20 1868   51]
 [  20 1887   73]
 [  20 1899   80]
 ...
 [ 436 2000  853]
 [ 436 2023 1001]
 [ 436 2042  917]]

Yet it now appears like this:

current transformed [[ -13.8649482   531.65490952  923.63328138]
 [ -13.89809115  580.32768386  977.55101852]
 [  -4.96404728  924.18309313  357.72738699]
 ...
 [ 446.53994182 1935.10420533 1131.60909333]
 [ 451.14326026 1943.64447165  637.97002858]
 [ 444.92508985 1965.00453948 1338.77464399]]