Add nD support to several geometric transform classes #5188
Conversation
|
Hello @jni! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:
Comment last updated at 2021-03-11 11:12:33 UTC |
|
The CI failure appears to be an unrelated hang. |
skimage/transform/_geometric.py
Outdated
| matrix = np.eye(dimensionality + 1) | ||
| else: | ||
| dimensionality = matrix.shape[0] - 1 | ||
| if matrix.shape != (dimensionality + 1, dimensionality + 1): |
There was a problem hiding this comment.
Shouldn't dimensionality be deduced from matrix.shape if matrix is not None?
I mean, when matrix is provided, we can assume that the user knows the dimensionality, we don't need to force him to give it as argument if dimensionality != 2...
I think dimensionality should be None by default...
There was a problem hiding this comment.
Well, sorry, I need one more coffee this morning 😅...
BTW, this situation can never happen since dimensionality is overwritten when matrix is not None.
There was a problem hiding this comment.
Waouh, I need a bucket of coffee...
There was a problem hiding this comment.
😂 I actually checked the diff cover when developing tests here so I'm fairly confident that all lines are reachable. 😃
skimage/transform/_geometric.py
Outdated
| nparam = v.size | ||
| # solve for d in: d * (d - 1) = nparam | ||
| d = (1 + np.sqrt(1 + 4 * nparam)) / 2 - 1 | ||
| dimensionality = int(d) |
There was a problem hiding this comment.
for nparam == 9, solution is, e.g., 3.999 — is int then the right call?
There was a problem hiding this comment.
Well, there's no such thing as 9 parameters for the affine... It's 6 (2D), 12 (3D), 20 (4D), ... And anyway you've made a mistake typing stuff in somewhere because I get 4.0 for n=20 and 3.88748 for 19?
But, it's totally true that one should never trust floating point arithmetic to be exact. I'll use np.round here.
skimage/transform/_geometric.py
Outdated
| v, (dimensionality, dimensionality + 1) | ||
| ) | ||
| matrix = np.concatenate( | ||
| (part_matrix, np.eye(dimensionality + 1)[-1:]), axis=0 |
There was a problem hiding this comment.
Easier to just pop the one into place?
There was a problem hiding this comment.
🤷 The alternative code is:
matrix = np.concatenate(
(part_matrix, np.zeros(dimensionality + 1), axis=0
)
matrix[-1, -1] = 1So, you get an extra line of code, for something maybe marginally clearer? I agree that what I'm after really is the vector commonly called e_i. But I don't think that exists in NumPy. =)
There was a problem hiding this comment.
It actually exists: what you are looking for is np.eye(1, dimensionality + 1, dimensionality) 😉
>>> np.eye(1, 4, 3)
array([[0., 0., 0., 1.]])| centered = points - centroid | ||
| rms = np.sqrt(np.sum(centered ** 2) / n) | ||
|
|
||
| if rms == 0: |
There was a problem hiding this comment.
I guess this point could come pretty close to zero by itself, so is it worth doing the explicit check, or should we just let the algorithm run its course?
There was a problem hiding this comment.
Or should we define a tolerance parameter for numerical stability?
if rms < tol:
...There was a problem hiding this comment.
tol is hard to specify, so I'd say "no, if the algorithm produces nans that pass through anyway" or "yes, if @jni says we need to force the output to have a certain context". Either way, there may be better ways to check than looking at the RMS---like looking at the condition number of the matrix being solved.
There was a problem hiding this comment.
OK, after reading more closely, I think that @jni wants to avoid these try... except... blocks
bingo! 😊 And those are only triggered when actual 0 is reached. I'll add a comment.
|
Ok @stefanv I think I've addressed all your questions! Thanks for the review! 🎉 |
| # We need the axes other than the rotation axis, in the right order: | ||
| # 0 -> (1, 2); 1 -> (0, 2); 2 -> (0, 1). | ||
| axes = sorted({0, 1, 2} - {axis}) | ||
| # We then embed the 2-axis rotation matrix into the full matrix. | ||
| # (1, 2) -> R[1:3:1, 1:3:1] = R2, (0, 2) -> R[0:3:2, 0:3:2] = R2, etc. | ||
| sl = slice(axes[0], axes[1] + 1, axes[1] - axes[0]) |
There was a problem hiding this comment.
Give me more magic
| # We need the axes other than the rotation axis, in the right order: | |
| # 0 -> (1, 2); 1 -> (0, 2); 2 -> (0, 1). | |
| axes = sorted({0, 1, 2} - {axis}) | |
| # We then embed the 2-axis rotation matrix into the full matrix. | |
| # (1, 2) -> R[1:3:1, 1:3:1] = R2, (0, 2) -> R[0:3:2, 0:3:2] = R2, etc. | |
| sl = slice(axes[0], axes[1] + 1, axes[1] - axes[0]) | |
| args = [None, None, None] | |
| args[-axis] = 1 + ( axis > 0) | |
| sl = slice(*args) |
😈
There was a problem hiding this comment.
That does look correct, but is pretty obscure (thus the devil emoji, I suppose!).
sorted axes could also be created with a generator expression:
# Get the (sorted) axes of the plane in which rotation occurs
ax1, ax2 = (ax for ax in range(3) if ax != axis)
# We then embed the 2-axis rotation matrix into the full matrix.
sl = slice(ax1, ax2 + 1, ax2 - ax1)There was a problem hiding this comment.
Imma keep my code + comments, thanks both. 😂
skimage/transform/_geometric.py
Outdated
| The Euler rotation matrix. | ||
| """ | ||
| if axes is None: | ||
| axes = list(range(3)) |
There was a problem hiding this comment.
| axes = list(range(3)) | |
| axes = tuple(range(3)) |
It doesn't look like the list is modified later, so may as well use a tuple
There was a problem hiding this comment.
There was a problem hiding this comment.
actually, we only iterate over axes later, so I'm just gonna leave it as range()!
skimage/transform/_geometric.py
Outdated
| Implemented only for 2D and 3D. For 3D, this is given in XZX Euler | ||
| angles. |
There was a problem hiding this comment.
XZX is not right here is it? It looks like we use axes in the order 0, 1, 2 in _euler_matrix which should be XYZ (assuming axis=0 means X?). XZX is apparently one of 12 possible conventions, but it didn't seem to be the one used here.
Another form for 3D that I find fairly intuitive is where the user specifies a unit vector corresponding to the [axis of rotation and then a single rotation angle about that vector. There is an explicit formula for that form.
There was a problem hiding this comment.
XZX is not right here is it?
Ah, yes, XZX is left over from an older version of the code. Will fix.
Another form for 3D that I find fairly intuitive is where the user specifies a unit vector corresponding to the axis of rotation and then a single rotation angle about that vector.
Could change to that if people think that's better...
There was a problem hiding this comment.
@grlee77 since the 3D Euler angles and the 3D vector + rotation are unambiguously distinguishable, I'm going to leave it as 3D Euler angles for now (fixed the docstring), and we can revisit in the future if someone clamors for alternate parameterizations.
This reverts commit 9630425. This change turned out to be incorrect because of the missing `norm_factor`.
|
@scikit-image/core I believe I've addressed all the review comments — please have another look! 🎉 |
|
Note, the optional deps failure is from #5266 and unrelated to this PR. |
|
Nice work. Thank you @jni! |
This is the second PR with more gradual changes from #3544. I've pulled out the
automatic conversion of 1D parameter vectors to transforms. Viewed in
isolation, that was definitely way too magical given our values. I've left the
implementation in as a private function for future reference. Eventually, I'd
like to add a class method,
from_parameter_vector, that would bewell-documented for each class, but I think that should come in a future PR.
I've done my best with the Euler rotation matrix code, but that would benefit
the most from review. Since we don't have a strict xyz/zyx axis order, it's
basically ambiguous whether XYZ or ZYX Euler angles are implemented (if it's
either of those 😂). CC @stefanv @grlee77 @matthew-brett
Description
Checklist
./doc/examples(new features only)./benchmarks, if your changes aren't covered by anexisting benchmark
For reviewers
later.
__init__.py.doc/release/release_dev.rst.