Add ellipsoid_coords function #5263
Conversation
|
Hello @ksugar! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:
Comment last updated at 2021-03-16 15:29:13 UTC |
based on the following review: scikit-image#5263 (comment) Co-authored-by: Juan Nunez-Iglesias <juan.nunez-iglesias@monash.edu>
|
Thanks for the update @ksugar! Now all that's missing is updating the tests to match 😬 😂, and as far as I'm concerned it's ready to go! 🎉 |
|
Hi @jni , thank you for your review! I am working on updating the docstrings and the tests. Could you wait a little until I update the code in that way? Thanks! |
|
@ksugar all good. I wonder though, whether a rotation matrix is too many degrees of freedom — indeed, a matrix may or may not be a rotation matrix. I'm wondering whether we should support both APIs: |
|
@jni You have a point! Then, how about the interface like this? , where 📝 Additional note: The motivation to support |
|
Hi @ksugar just a question about the API: I'm not sure why we need two functions here, |
|
Hi @emmanuelle |
|
@emmanuelle @jni Please give me advice or comments on whether it would be better to combine |
|
@ksugar you're absolutely right that ellipsoid was the wrong API for draw, and I don't think we should harmonize the APIs. I think medium term we can deprecate ellipsoid and replace it with ellipsoid coords, but that deprecation can happen after this PR. @emmanuelle do you agree?
My suggestion is that we leave the Other than that, |
|
@jni thank you for your feedback and suggestions! I see your point that it is better to start from the minimal implementation. I will leave the I found that we need to specify the order of rotations (e.g. |
|
Wow, lots of updates, thanks @ksugar! 😅 Some further comments:
|
|
I 👍 @jni suggested changes, I would also make |
|
Hi @jni thank you for your feedback! I would like to confirm about the will generate works perfectly in my case, where the transpose operation I would also like to use the term |
|
Hahaha "your function is amazing! btw, it is also wrong." 😂 I think you are almost certainly right @ksugar! 🤦 Would be great if @grlee77, @stefanv or @rfezzani could check the above, as I'm a total newcomer to 3D rotations and they make my head hurt. 😕 But by my reading, we did get it wrong in #5188, even after much discussion. 😂 Thanks @ksugar! |
|
@jni actually your code is great 👍 I'm learning a lot of techniques from your code! If this change is approved by the members and to |
|
kind ping to @stefanv @grlee77 @rfezzani, maybe @JDWarner, if someone could check our math for 3D Euler rotations it would be much appreciated! See discussion starting at #5263 (comment) |
|
At first, I agreed with @ksugar, rotations should applied from left to right. But to convince myself, I wanted to compare >>> import copy
>>> import open3d as o3d
>>> from skimage.transform._geometric import _euler_rotation_matrix
>>> angles = (np.pi/2,0,np.pi/4)
>>> mesh = o3d.geometry.TriangleMesh.create_coordinate_frame()
>>> mesh_r_o3d = copy.deepcopy(mesh).translate((2,0,0))
>>> mesh_r_skimage = copy.deepcopy(mesh).translate((4,0,0))
>>> mesh_r_o3d.rotate(mesh.get_rotation_matrix_from_xyz(angles))
TriangleMesh with 1134 points and 2240 triangles.
>>> mesh_r_skimage.rotate(_euler_rotation_matrix(angles))
TriangleMesh with 1134 points and 2240 triangles.
>>> o3d.visualization.draw_geometries([mesh, mesh_r_o3d, mesh_r_skimage])
>>> np.array_equmesh.get_rotation_matrix_from_xyz(angles)
np.array_equal( np.array_equiv(
>>> np.array_equal(mesh.get_rotation_matrix_from_xyz(angles), _euler_rotation_matrix(angles))
TrueThe results are equal! |
|
@rfezzani @jni it is great to compare the output with Open3D! I checked the source code of Open3D and its implementations are exactly same as what we can find in the I would like to confirm the assumptions in the implementations.
If On the other hand, I was assuming that the functions rotate standard basis vectors and the ordering is About the order of matrices in calculation, as far as I understand, the following explains it. If we define |
|
I don't have the background knowledge to comment on correctness. @jni are you able to do one more round of review so we can get this in? I wonder if one of our resident scipy rotation experts—@pmla, @nmayorov, @adbugger, or @evbernardes—can help out? |
|
@evbernardes regarding the following point, what I thought was drawn in the attached figure.
If we rotate the axes with a rotation matrix
|
|
@ksugar if understood correctly, it depends how you are defining the rotation matrix. A rotation matrix is just a representation of a rotation linear transformation, that could be define either from the inertial frame to the body frame, or from the body frame to the inertial frame. I think by putting an inversion inside of the implementation could add some confusion, and I'd argue this is one of the reasons while just taking a matrix directly (or a Scipy Rotation instance) is a better solution. Edit.: Of course, it depends if a certain standard is assumed for this application, but I can't comment on that since I'm not familiar with it! |
|
I'm more used to defining the rotations from the body frame to the inertial frame. This means that, for example, defining a normal vector to some body along the body z axis, it would always look like (0 , 0, 1) in the body frame, and applying the rotation to it would express how you'd see it in the inertial frame. |
|
I forgot to mention: if the rotation matrix as the linear transformation that transforms one vector to the rotated version of the same vector, a unit vector that defines an axis is still just a vector, so there's nothing special about it, I think. And to rotate a point, depends also on how you define the point, but points can also be described by displacement vectors. Hope that makes sense with your application! |
draw3d.ellipsoid_coords is a function analogous to draw.ellipse. This function returns the voxel coordinates of an ellipsoid, while the existing function draw3d.ellipsoid returns a 3d ndarray that renders an ellipsoid.
based on the following review: scikit-image#5263 (comment) Co-authored-by: Juan Nunez-Iglesias <juan.nunez-iglesias@monash.edu>
- properly formulate rotations based on rotation_angles, rotation_order - rotation_matrix support - fix and add test cases
- remove rotation_matrix from arguments of ellipsoid_coords
- update to use skimage.transform._geometric._euler_rotation() to generate a rotation matrix - update docstrings and tests - use .T instead of np.linalg.inv to compute a inverse of a rotation matrix - use R to represent a rotation matrix
This commit is a part of the last commit. - update to use skimage.transform._geometric._euler_rotation() to generate a rotation matrix - update docstrings - use .T instead of np.linalg.inv to compute a inverse of a rotation matrix - use R to represent a rotation matrix
stefanv
force-pushed
the
ellipsoid-coords
branch
from
529f32c
to
dba5fb2
Compare
|
@ksugar Are you happy with where this is at, given @evbernardes's comments above? I think we are ready to merge, and can get this into the 0.21 release. |
- use scipy.spatial.transform.Rotation.from_eular to create a rotation matrix in draw3d.ellipsoid_coords()
|
@stefanv thank you for the updates and I'm glad to hear that this PR is going to be merged soon! |
|
This could make a nice example for the gallery, if you're up for it. Meanwhile, I'll push a few touch-ups and get this merged. |
|
Ran into some more questions; also @jni may have thoughts:
d_lim, r_lim, c_lim = np.ogrid[:float(bounding_shape[0]),
:float(bounding_shape[1]),
:float(bounding_shape[2])]
d_org, r_org, c_org = scaled_center - upper_left_bottom
d_rad, r_rad, c_rad = axis_lengths
conversion_matrix = R.T @ np.diag(spacing)
d, r, c = (d_lim - d_org), (r_lim - r_org), (c_lim - c_org)
distances = (
((d * conversion_matrix[0, 0]
+ r * conversion_matrix[0, 1]
+ c * conversion_matrix[0, 2]) / d_rad) ** 2 +
((d * conversion_matrix[1, 0]
+ r * conversion_matrix[1, 1]
+ c * conversion_matrix[1, 2]) / r_rad) ** 2 +
((d * conversion_matrix[2, 0]
+ r * conversion_matrix[2, 1]
+ c * conversion_matrix[2, 2]) / c_rad) ** 2
)Ideally, we would not be repeating the code for three axes. So, while we're close on this PR, I'm going to take this off the 0.21 milestone, and we can get it into 0.22. We're on an increased release cadence now, so shouldn't take too long. |
|
Happy I was able to help! Here are some extra thoughts:
I think there are two advantages of using axis by their names
If the focus on Scikits is to be extensions of SciPy, I'd argue that consistency with SciPy should be a priority over consistency with NumPy, don't you think? |
|
|
||
| # Generate a rotation matrix. The order of the elements needs to be | ||
| # reversed so that it follows the (z, y, x) order. | ||
| R = Rotation.from_euler(seq, angles).as_matrix()[::-1, ::-1] |
There was a problem hiding this comment.
What is this inversion doing exactly? And why is axes_str defined as ZYX instead of XYZ?
Are there two definitions just inverting one another?
There was a problem hiding this comment.
Typically, axis 0 corresponds not to X but to Z. Think, e.g., of an image represented as an array: (M, N) means (Y, X), upside down, sort of. You see the confusion :)
There was a problem hiding this comment.
I see, this is because of the application for images then!
| :float(bounding_shape[2])] | ||
| d_org, r_org, c_org = scaled_center - upper_left_bottom | ||
| d_rad, r_rad, c_rad = axis_lengths | ||
| conversion_matrix = R.T @ np.diag(spacing) |
There was a problem hiding this comment.
And why is the inverse of the rotation used here?
There was a problem hiding this comment.
I woke up last night, realizing that I had lied to you. It gets more complicated: scikit-image, in a regrettable but somewhat rational early API decision, chose to adopt the XYZ coordinate system for transformations. So that's layered on top of everything. We'll fix that in v2.
There was a problem hiding this comment.
I woke up last night, realizing that I had lied to you.
Good to know this doesn't only happen to me!
There was a problem hiding this comment.
Here's an idea, for a future implementation.
What about a new scikit.Rotation class that inherits from scipy.spatial.transform.Rotation class, but only changes the apply method in order to invert which axes to use (XYZ or ZXY), and as_euler/from_euler to use your preferred standard of 0,1,2 instead of X,Y,Z?
This shouldn't be too time-consuming to implement with inheritance and would allow everyone to use whatever they prefer using (matrices, quaternions or angles).
| f'len(axis_lengths) should be 3 but got {len(axis_lengths)}') | ||
| axis_lengths = np.array(axis_lengths) | ||
|
|
||
| if angles is None: |
There was a problem hiding this comment.
Btw, you don't have to use rotation matrices. You can just define the Rotation object directly and use Rotation.apply() which is a lot more efficient.
I think SciPy is making a mistake here. When working with arrays, we understand what axes 0, 1, and 2 are; but after two decades if doing that, if you ask me which ones represent X, Y, and Z, I wouldn't know. It depends on how you view your array? |
I see, personally I'd say it's best to have it consistently, since they are most likely being used together. But ofc, this is up to you |
|
I still think it's best to use rotation matrices directly (or the rotation object) as an input though, since it introduces a possible very inefficient step of first decomposing the rotation into ZYX angles and then converting it to rotation matrices again |
Description
This PR adds the
draw3d.ellipsoid_coordsfunction that is analogous to thedraw.ellipsefunction.This function returns the voxel coordinates of an ellipsoid, while the existing function
draw3d.ellipsoidreturns the 3d ndarray that renders an ellipsoid.This function accepts the rotations in each axis (i.e.
rot_x,rot_yandrot_z).It also supports the anisotropic resolution by introducing the
spacingparameter.A variation of this function is implemented and used in the following preprint.
Original implementation used in the preprint:
https://github.com/elephant-track/elephant-server/blob/v0.1.0/elephant-core/elephant/util/ellipsoid.py
Checklist
Docstrings for all functions
Gallery example in
./doc/examples(new features only)I am not sure if I should add the example.
Benchmark in
./benchmarks, if your changes aren't covered by anexisting benchmark
not applicable
Unit tests
Clean style in the spirit of PEP8
For reviewers
later.
__init__.py.doc/release/release_dev.rst.