Affine transform: Shearing not implemented

[haesleinhuepf]

In the current implementation of affine transforms, shearing has not been implemented. I don't use shearing that's why I'm not sure how the API should look like. In CLIJ2 affine transforms, we used that terminology:

shearXY=[factor]: shearing along X-axis in XY plane according to given factor
shearXZ=[factor]: shearing along X-axis in XZ plane according to given factor
shearYX=[factor]: shearing along Y-axis in XY plane according to given factor
shearYZ=[factor]: shearing along Y-axis in YZ plane according to given factor
shearZX=[factor]: shearing along Z-axis in XZ plane according to given factor
shearZY=[factor]: shearing along Z-axis in YZ plane according to given factor

I'm happy to implement this here again, but I'd need someone who tests it and provides feedback.

[tlambert03]

Could you use an axes argument (tuple of 2 ints), as with done in the scipy.ndimage API? https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.rotate.html#scipy.ndimage.rotate

[haesleinhuepf]

Hey @tlambert03 ,

this link points to rotation, not shearing. And I can't find shearing in scipy...

Anyway, I think I'll take a closer look on affine transforms in skimage. It might make sense to make the API compatible with this one.

[haesleinhuepf]

... and I just realized that skimage only supports 2D :-(

[tlambert03]

Correct scipy.ndimage.interpolate doesn't offer a shear convenience function (just zoom, rotate, and shift ... so I pointed to rotate as an example). But the concept is the same: axes is a two-tuple of ints that defines the shear plane, (and ordering of provided axis matters to determine the shearing axis). factor then would be a second argument that determines the magnitude.

... and yeah, skimage only supports 2D.

I'm pretty sure I've mentioned this before, but you can have a look at my implementation at ndtiffs. Though it's a subset of functionality, I tried pretty hard to match the scipy.ndimage API there for what I did implement, and made tests to make sure that the two APIs yielded similar results. (i.e. that repo is not just copy-paste from gputools... I changed a lot to make it more like other python apis)

[tlambert03]

you can see i stopped short of provided the shearing axis as an input argument, and just hard-coded it here in the get_deskew_func function... but the (2, 0) there (in mat[2, 0] = -deskewFactor) is where you'd put the axes argument to build a matrix to pass to the full affine function.

[tlambert03]

so basically ... my vote would be to make this all a drop-in replacement for the scipy.ndimage.interpolate API as much as possible. And if you'd like to offer a shear convenience function (which they do not) just mimic the rotate function with regards to axis declaration. @jni have talked about wanting a drop-in replacement like that, and ndtiffs was starting to offer that (see the try: import here where it simply gracefully falls back to scipy.ndimage if pyopencl isn't installed... without changing the user API)

[haesleinhuepf]

Just for the record, I'm adding the scipy affine_transform signature here:

scipy.ndimage.affine_transform(
    input, 
    matrix, 
    offset=0.0, 
    output_shape=None, 
    output=None,
    order=3, 
    mode='constant', 
    cval=0.0, 
    prefilter=True
)

and the clesperanto version

pyclesperanto_prototype.affine_transform(
    source : Image, 
    destination : Image = None, 
    transform : Union[np.ndarray, AffineTransform3D, AffineTransform] = None, 
    linear_interpolation : bool = False
)

I'll need to think a bit about how to achieve a "drop in replacement" without breaking the interoperability with the Java side. It's tricky ;-)

[tlambert03]

I'll need to think a bit about how to achieve a "drop in replacement" without breaking the interoperability with the Java side. It's tricky ;-)

It's tricki_er_ for sure... but doable. and as I mentioned over in #91 ... I think this really gets to the heart of what kind of bridge you want to build :) One is "learn the clij way and you can use it anywhere" (which has its merits!), the other is "do it the natural way in your ecosystem ... and clij will take care of building the bridge"

[haesleinhuepf]

The good thing is: We're changing API-bits anyway when switching from clij to clesperanto, also on Java side. Thus, we can change the API. I just would like to make sure that these functions are no exceptional cases. The should fit in the general API style of clesperanto...

In other words: All we change here, should also be possible on Java side. Otherwise it'll be a one-way bridge ;-)

[pr4deepr]

Hi @haesleinhuepf
I am happy to have a look and provide feedback.
I quite like the idea of having the terminology similar to CLIJ, where you define the shear factor. It was extremely convenient and approachable to use strings and enter values to perform shearing and other affine transformations, compared to specifying transformation matrices. Of course, having a few examples and some documentation would be handy.

For context, I was using CLIJ to visualise the affine transformations using the CLIJ-assistant in FIJI and then export the same code to use in Python. I guess the syntax for shearing you want to adopt will depend on the audience you are targeting as well.

Cheers
Pradeep

[pr4deepr]

Just leaving this here to continue our discussion for implementing shearing..
I haven't tested it much..

In regard to the shearing, the matrix will be of the form:

I think we can implement shearing in x and y direction by calculating deskew factor like you did:

deskew_factor = 1.0 / math.tan(rotation_angle * math.pi / 180)

For shear in z along Y direction
we change h(yz) I think.

Similarly for shear in z along X, we change h(xz).

I think shear operation may look something like this in clesperanto:

    def shear_z(self,rotation_angle_x: float = 0, rotation_angle_y: float = 0 ):
   
        import math
        try:
            deskew_factor_x = 1.0 / math.tan(rotation_angle_x * math.pi / 180)
        except ZeroDivisionError:
            deskew_factor_x = 0
        
        try:
            deskew_factor_y = 1.0 / math.tan(rotation_angle_y * math.pi / 180)
        except ZeroDivisionError:
            deskew_factor_y = 0


        # shearing
        self._concatenate(np.asarray([
            [1, 0, deskew_factor_x, 0],
            [0, 1, deskew_factor_y, 0],
            [0, 0 , 1, 0],
            [0, 0, 0, 1],
        ]))
        return

This is only for shear in z along the x and y direction. These are the two major ones..

Won't get time to test this today though. Will get onto this soon..

Source for matrix: https://people.cs.clemson.edu/~dhouse/courses/401/notes/affines-matrices.pdf

[haesleinhuepf]

Awesome, thanks @pr4deepr! You're shearing/deskewing data more often than I do, that's why I'm super excited that you proposed this API. I looks very good to me. If it works for you, feel free to submit a PR. If we know how to make it useful/handy for one plane, I'm also happy to implement it for the other planes. Let me know how it goes and if you need anything! And of course: no hurry :-)

[haesleinhuepf]

Jeey, this has been solved with the fantastic support by @pr4deepr in #156 and #157. 🥳