registration interpolation

[orena1]

Hi @glichtner
I'm trying to understand the difference between using skimage wrap and pystackreg.
From the documentation of skimage - https://scikit-image.org/docs/stable/api/skimage.transform.html#skimage.transform.warp I see that the default interpolation is Bi-linear, what is the interpolation used for pystackreg? (Assuming Affine transformation).

also, correct me if I'm wrong, but can you directly use the tmat in skimage wrap function?

for i in range(tmats.shape[0]):
    #tform = tf.AffineTransform(matrix=tmats[i, :, :])
    out[i, :, :] = tf.warp(img1[i, :, :], tmats[i, :, :]) # tform)

Thanks!

[glichtner]

Hi @orena1,

pystackreg uses cubic spline interpolation for transforming images.

Regarding the tf.warp function: I think you are right and you can use the tmat directly, but I haven't tried it yet - have you?

Best,

[orena1]

Thanks @glichtner ,

pystackreg uses cubic spline interpolation for transforming images.

It does not seem to give me identical images...

import math
import numpy as np
import matplotlib.pyplot as plt

from skimage import data
from skimage import transform
from skimage import img_as_float
from pystackreg import StackReg
sr = StackReg(StackReg.AFFINE)


tmat = np.array([[  1.02,   0.  , -17.29],
       [  0.  ,   0.97,  -0.49],
       [  0.  ,   0.  ,   1.  ]])


img = img_as_float(data.chelsea())
img = img[:,:,0] # get only first channel
img_sk = transform.warp(img, tmat,order=3) ## order 3 is cubic https://scikit-image.org/docs/dev/api/skimage.transform.html#skimage.transform.warp
img_tr = sr.transform(img,tmat=tmat)

np.all((tf_img-img_tr)==0)
Out[1]: False

plt.imshow(tf_img - img_tr,cmap='gray')

---

Regarding the tf.warp function: I think you are right and you can use the tmat directly, but I haven't tried it yet - have you?

It seems to be identical

import math
import numpy as np
import matplotlib.pyplot as plt

from skimage import data
from skimage import transform
from skimage import img_as_float
from pystackreg import StackReg
sr = StackReg(StackReg.AFFINE)
tmat = np.array([[  1.02,   0.  , -17.29],
       [  0.  ,   0.97,  -0.49],
       [  0.  ,   0.  ,   1.  ]])

img = img_as_float(data.chelsea())
img = img[:,:,0] # get only first channel

img_sk0 = transform.warp(img, tmat,order=3)

tform = transform.AffineTransform(matrix=tmat)
img_sk1 = transform.warp(img, tform,order=3)
np.all((img_sk0-img_sk1)==0)
Out[1]: True

[glichtner]

Hi @orena1,

Thanks for the detailed testing and the pull request.
I don't know what the exact reason for the difference between skimage's implementation of the cubic spline interpolation (I think they use spline interpolation, but it's not clear to me from the docs) and the one from turboreg is, but any difference in data types (e.g. float precision), rounding operations, order of calculations and whatnot can lead to numerical differences in the result and therefore I am not worried by the difference and neither do I think the difference is particularly important.

[orena1]

Thanks

[patquem]

Hi,
I was faced to the same problem/question as @orena1
When regarding the range of my registered image, I was surprised to meet out of range data for a simple translation transformation.
I am not sure that considering a 3rd order interpolation is a good idea for the transformation.
I would recommand a 1rst order, with the possibility given to the user to increase the order of interpolation (at his own risk).
As for @orena1 , when comparing the results returned by pystackreg with a 3rd order interpolation from skimage.transform.warp, the results differ.
The out of range results returned by pystackreg interpolation is on my opinion not a good feature (by default) :(
Patrick

import numpy as np
from scipy import ndimage as ndi
from skimage.data import binary_blobs
from skimage.transform import warp, AffineTransform
from pystackreg import StackReg

ref = binary_blobs(300,
                   blob_size_fraction=0.05,
                   volume_fraction=0.7,
                   seed=0).astype(float)

shift = [6.2, 3.8]
mov = ndi.shift(ref, shift, order=1)

# cropping
ref = ref[20:-20, 20:-20]
mov = mov[20:-20, 20:-20]

# Translational transformation
sr = StackReg(StackReg.TRANSLATION)
out_tra1 = sr.register_transform(ref, mov)

tform = AffineTransform(matrix=sr.get_matrix())
out_tra2 = warp(ref, tform, mode='constant', cval=0.)
out_tra3 = warp(ref, tform, order=3)

print(np.max(out_tra1)) # -> 1.0594804286956787
print(np.max(out_tra2)) # -> 1.0
print(np.max(out_tra3)) # -> 1.0

[glichtner]

Dear @patquem,

Thank you for your comment and the example code. I agree that out of range results are not desirable, but for the time being I had decided to reproduce exactly the results from the ImageJ plugins TurboReg/StackReg, which show the same behavior. I will consider changing that behavior in the future and possibly rely on the skimage package to perform the transformation instead of using the original code from TurboReg/StackReg.

Best,
Gregor

[patquem]

ok @glichtner.
I understand your reasons.
As long as you haven't done the 'future' changes, maybe it would be good to advice the user in the documentation about what he can be faced and how he can go back to a lower and user defined interpolation thanks to skimage warp method (as suggested below), wouldn't it be ?
Thanks,
Patrick

from skimage.transform import AffineTransform, warp 
tform = AffineTransform(matrix=sr.get_matrix())
registered_img = warp(img, tform, order=1, mode='constant', cval=0.)