raphael-group / paste

Probabilistic Alignment of Spatial Transcriptomics Experiments
https://paste-bio.readthedocs.io/en/latest/
BSD 3-Clause "New" or "Revised" License
77 stars 18 forks source link

Extracting the rotation matrix from `stack_slices_pairwise` #21

Closed leoguignard closed 2 years ago

leoguignard commented 2 years ago

Hello, first thanks a lot for PASTE!

While using your algorithm to align two slices, it can be interesting to filter out beads/cells that do not have enough expression (as you do for some of your examples).

That being said, it can also be useful to apply the transformation computed to the original slice (with all the beads/cells). One way to do so is to get the computed transformations from generalized_procrustes_analysis. At the time being, stack_slices_pairwise returns the translation and the rotation angle computed from the rotation matrix. I personally find that the rotation angle can be ambiguous (for example not knowing its origin or if there were a flip in the transformation). This is why I suggest that stack_slices_pairwise could also return the rotation matrix instead of only the angle. I have implemented such a feature in the fork that you can find there. The change suggested in the fork should not modify the original behavior, it just allows the user to have access to the rotation matrix instead of the angle if wanted.

Then, using the output rotation together with the translation, one can then apply it to any other set of points for example using the following function:

def apply_trsf(M, trans, points):
    """
    Apply a rotation from a 2x2 rotation matrix `M` together with
    a translation from a translation vector of length 2 `trans` to a list of
    `points`

    Args:
        M (nd.array): a 2x2 rotation matrix
        trans (nd.array): a translation vector of length 2
        points (nd.array): a nx2 array of `n` points 2D positions

    Returns:
        (nd.array) a nx2 matrix of the `n` points transformed
    """
    trsf = np.identity(3)
    trsf[:-1, :-1] = M
    tr = np.identity(3)
    tr[:-1, -1] = -trans
    trsf = trsf @ tr

    flo = points.T
    flo_pad = np.pad(flo, ((0, 1), (0, 0)), constant_values=1)
    return ((trsf @ flo_pad)[:-1]).T

Note that stack_slices_pairwise could also directly return the homogeneous transformation (trsf in the previous function) so that one can directly apply it to a set of points also in homogeneous coordinates.

I would be happy to open a PR for the previously mentioned change if you would be interested in it. Moreover, if you think the whole 4x4 transformation matrix would be better than just the rotation one, I would be happy to propose another PR for that.

Have a nice day, Cheers, Léo

mrland99 commented 2 years ago

Hey Léo,

Sorry for the delayed response. I think returning the rotation matrix is a great idea! Definitely make a PR for that.

And also please make another PR for your apply_trsf() function, maybe underhelper.py?

I don't know if the whole 4x4 transformation is necessary (unless you think there are more use cases than just the rotation)-- I think just the rotation is good, especially since the change is pretty easy and straightforward.

But yeah, thanks for doing this and let me know if you have any other ideas but your current proposed implementation sounds good to me.

Cheers, Max

leoguignard commented 2 years ago

Hey @mrland99,

Thanks for the answer, of course no problem for the delay :)

How do you feel about having only one PR for both the return of the rotation matrix and the apply_trsf function un helper.py? I'd rather do only one PR if it's fine with you.

Otherwise, if you'd rather keep them split, I definitely could understand too :)

Cheers, Léo

mrland99 commented 2 years ago

One PR works if it's easier for you. Thanks Léo.

mrland99 commented 2 years ago

Merged and published. Thanks again for the feature-- I made sure to give you a shoutout in the release notes for the latest version :)