Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixed a bug in __mul__ #24

Closed
wants to merge 1 commit into from
Closed

Conversation

ToddSmall
Copy link

I noticed that composition of transformations doesn't work, at least some of the time. Here's an example:

In [1]: import affine            

In [2]: affine.__version__
Out[2]: '1.2.0'

In [3]: from affine import Affine

In [4]: Affine.rotation(90.) * (Affine.translation(-10., -5.) * (12, 5))      
Out[4]: (0.0, 2.0)

In [5]: (Affine.rotation(90.) * Affine.translation(-10., -5.)) * (12, 5)
Out[5]: (-10.0, 22.0)

I think the problem is due to sign flips in the __mul__ method. I've corrected the sign flips where necessary and added some additional test coverage. All tests pass.

However, it's possible that I'm simply not understanding some nuance of raster coordinates versus "standard" Cartesian coordinates.

(ca, sa, 0.0,
-sa, ca, 0.0,
(ca, -sa, 0.0,
sa, ca, 0.0,
0.0, 0.0, 1.0))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems that this modification simply switches from clockwise convention to counter-clockwise. See https://en.wikipedia.org/wiki/Transformation_matrix#Rotation

@mwtoews
Copy link
Contributor

mwtoews commented May 2, 2016

I don't see any problems with the examples. What were the expected results from each?

@ToddSmall
Copy link
Author

The key point is that the two examples should have the same result. The first result, (0., 2.), is correct.

@@ -385,7 +385,7 @@ def __mul__(self, other):
vx, vy = other
except Exception:
return NotImplemented
return (vx * sa + vy * sd + sc, vx * sb + vy * se + sf)
return (vx * sa + vy * sb + sc, vx * sd + vy * se + sf)
Copy link
Contributor

@mwtoews mwtoews May 3, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks like a bug with the correct fix. The 3x3 matrix is [sa, sb, sc; sd, se, sf; 0, 0, 1] and the 3x1 matrix is [vx; vy; 1]. Can anyone else confirm?

@coveralls
Copy link

coveralls commented May 3, 2016

Coverage Status

Coverage increased (+0.6%) to 98.052% when pulling 1d9c0d4 on ToddSmall:master into fca54b9 on sgillies:master.

@mwtoews
Copy link
Contributor

mwtoews commented May 3, 2016

I've given it a better look, and there might be a bug with the 3x3 dot 3x1 part of mul. Both results would then give (0.0, -2.0). Here is some proof with NumPy:

import numpy as np
# R = rotation, T = translation, P = point
R = np.array(Affine.rotation(90.)).reshape((3, 3))
T = np.array(Affine.translation(-10., -5.)).reshape((3, 3))
P = np.array([12, 5, 1])

# Line 4
np.dot(R, np.dot(T, P))  # [ 0. -2.  1.]

# Line 5
np.dot(np.dot(R, T), P)  # [ 0. -2.  1.]

(Note, the 3rd item for points is always 1)

The result (0, -2) is expected with this thought exercise:

  • point (12, 5) is translated by (-10, -5) to make (2, 0)
  • point (2, 0) is then rotated clockwise 90 to (0, -2)

However, I'm not sure about the other changes to rotation, since they change the chirality of the module, which is a bit invasive. However, I don't see any docstrings for rotate to indicate if the intended convention is either clockwise or anti-clockwise.

@ToddSmall
Copy link
Author

Thanks @mwtoews for looking at this PR carefully. I apologize for not including more explanation (e.g., an example using NumPy like the one you constructed).

The changes to the sign of rotation are, as you say, invasive, but they are required (I think -- I need to double check) once you correct for the error in __mul__ for a +45 deg rotation of (1, 1) to yield (0, sqrt(2)), as shown in the readme file.

@mwtoews
Copy link
Contributor

mwtoews commented May 3, 2016

It seems that up till now it was using a anti-clockwise convention, which I can see with playing with a raster in QGIS. There are different conventions used used in mapping, see [Bearing](https://en.wikipedia.org/wiki/Bearing_%28navigation%29#Bearing measurement). An alternative to streamline this convention is to use a -ve rotation:

# using the version that fixes only the 3x3 dot 3x1 bug
Affine.rotation(-45.0) * (1.0, 1.0)  # (-3.3306690738754696e-16, 1.414213562373095)

which is effectively (0, sqrt(2)). I wonder if something like this should be used internally:

    @classmethod
    def rotation(cls, angle, pivot=None, clockwise=False):
...
        if not clockwise:
            angle *= -1.0

@perrygeo
Copy link

perrygeo commented May 3, 2016

@ToddSmall I think you've identified a bug. It seems to only affect transforming point coordinates from rotated affines (e.g. rotated_affine * (col, row) is wrong) and your change on line 388 looks like the fix. The fact that the calls were not associative (composable) was a symptom.

So I'd propose the following:

  • Remove the changes to rotation direction, agree with @mwtoews that it's somewhat unnecessary and a breaking change.
  • open a new ticket for an optional clockwise kwarg if you want to pursue that further
  • keep the test_composable but maybe rename to test_associative to stick with the matrix math terminology
  • add an additional case to directly test forward and backward coordinate transforms on a rotated affine
  • fix the README; I'm fairly certain that the Affine.rotation(45.0) * (1.0, 1.0) result is incorrect and should be changed to (1.414213562373095, 1.1102230246251565e-16). You can run py.test README.rst to run the doctests directly.

Thanks for digging into this! Affine's use being primarily with geospatial rasters, we don't often deal with rotation parameters so it's easy to see how that might've slipped though the cracks.

@sgillies
Copy link
Member

sgillies commented May 3, 2016

@ToddSmall thanks for bringing this! I'm with @perrygeo above on next steps and happy to schedule a new release ASAP.

@ToddSmall
Copy link
Author

Thanks for the feedback @perrygeo and @sgillies. I'll post a new patch set this morning (PDT).

@mwtoews
Copy link
Contributor

mwtoews commented May 3, 2016

@sgillies @perrygeo it wasn't clear what a default rotation should be. Keep it at clockwise=False (current behaviour) or change it to clockwise=True (and modify the README.rst example)?

ToddSmall added a commit to ToddSmall/affine that referenced this pull request May 3, 2016
in rasterio#24 (comment)

o I've changed the rotation direction back to what was there before.
o I've renamed test_composable to test_associate.
o I've added an addition roundtrip test.
o I've fixed the rotation example in README.rst.
@coveralls
Copy link

coveralls commented May 3, 2016

Coverage Status

Coverage increased (+0.6%) to 98.052% when pulling 5820e02 on ToddSmall:master into fca54b9 on sgillies:master.

@ToddSmall
Copy link
Author

@sgillies Would you prefer that I squash these two commits into one?

@sgillies
Copy link
Member

sgillies commented May 4, 2016

@ToddSmall @mwtoews let's keep the original angle sign, which was positive == counter-clockwise, and has been reversed in this PR.

@ToddSmall squashing would be great, thanks!

A couple of the affine parameters in the affine * vector path
in __mul__ were swapped. I also added some additional test
coverage and fixed the rotation example in README.rst.

I have *not* changed the rotation direction convention.
@coveralls
Copy link

coveralls commented May 4, 2016

Coverage Status

Coverage increased (+0.6%) to 98.052% when pulling 47cf6f0 on ToddSmall:master into fca54b9 on sgillies:master.

@ToddSmall ToddSmall changed the title Fixed a bug (?) in __mul__ Fixed a bug in __mul__ May 4, 2016
@sgillies
Copy link
Member

sgillies commented May 6, 2016

@ToddSmall this PR has changed the sense of the rotation angle and that still needs to be reverted.

To demonstrate and guard against future regressions (sign changes being pernicious), I've added this new test to the master branch: https://github.com/sgillies/affine/blob/master/affine/tests/test_rotation.py.

This test fails with a checkout of your branch:

$ python -m pytest affine/tests
============================================================== test session starts ==============================================================
platform darwin -- Python 3.5.1, pytest-2.9.1, py-1.4.31, pluggy-0.3.1
rootdir: /Users/sean/code/affine, inifile:
plugins: cov-2.2.1
collected 52 items

affine/tests/test_pickle.py ..
affine/tests/test_rotation.py F
affine/tests/test_transform.py .................................................

=================================================================== FAILURES ====================================================================
______________________________________________________________ test_rotation_angle ______________________________________________________________

    def test_rotation_angle():
        """A positive angle rotates a vector counter clockwise

        (1.0, 0.0):

            |
            |
            |
            |
            0---------*

        Affine.rotation(45.0) * (1.0, 0.0) == (0.707..., 0.707...)

            |
            |      *
            |
            |
            0----------
        """
        x, y = Affine.rotation(45.0) * (1.0, 0.0)
        assert round(x, 14) == round(sqrt(2.0) / 2.0, 14)
>       assert round(y, 14) == round(sqrt(2.0) / 2.0, 14)
E       assert -0.70710678118655 == 0.70710678118655
E        +  where -0.70710678118655 = round(-0.7071067811865475, 14)
E        +  and   0.70710678118655 = round((1.4142135623730951 / 2.0), 14)
E        +    where 1.4142135623730951 = sqrt(2.0)

affine/tests/test_rotation.py:27: AssertionError
====================================================== 1 failed, 51 passed in 0.15 seconds ======================================================

With your code, the unit vector has been rotated clockwise, to below the x axis.

@ToddSmall
Copy link
Author

ToddSmall commented May 6, 2016

@sgillies Sorry for the grief. I am confused about two things:

  1. Are you intending to use a right-handed or left-handed coordinate system? In your comment above, you've drawn the more widely used right-handed system, not a left-handed system as is often used for 2D computer graphics and image rasters.
  2. In @perrygeo's comment above, his 45 degree rotation rotates the vector clockwise in a right-handed system and counter-clockwise in a left-handed coordinate system.

Addendum:

In fact, my version does rotate counter-clockwise for a left-handed coordinate system.

@perrygeo
Copy link

perrygeo commented May 6, 2016

To clarify my logic here:

The intention of rotation has always been clockwise

>>> am = np.array(Affine.rotation(90.0)).reshape(3,3)
>>> point = np.array([10,0,1])
>>> np.matmul(am, point)
array([  0., -10.,   1.])

Point (10,0) becomes (0, -10): clockwise rotation.

The __mult__ operation for rotated_affine * (col, row) behaved counter-clockwise

>>> # affine 1.3.0
>>> Affine.rotation(90.0) * (10, 0)
(0.0, 10.0)

So one of the behaviors has to change if the math is to be consistent and associative.

Fixing the __mult__ is mathematically correct (It makes __mult__ consistent with manually multiplying the matrix produced by rotation). So my conclusion was that line 388 was a bug and all previous behavior of the mult operator wrt rotated affines should be assumed to be buggy as well. IOW the counter-clockwise rotation order is an accidental result of that bug. The order implied by the code seems to be clockwise since that is the behavior of the rotation matrix.

I could see the case for the opposite approach: Assume the mult operator always produced correct results (by accident) but the rotation matrix has a bug and should be corrected to counter-clockwise. This is effectively the approach that @ToddSmall took initially and may be better (less impact to users who rely on the rotation of the current mult operator).

In either case, the direction just needs to be made consistent to resolve this bug. The question is where do we break backwards compatibility?

@sgillies
Copy link
Member

sgillies commented May 14, 2016

We're disagreeing on which one of the two rotation conventions described in http://mathworld.wolfram.com/RotationMatrix.html we should follow. SVG, for example, follows the second convention: a rotation(90) appears to rotate an object clockwise, but in fact what it does is rotate the user coordinate system containing the object 90 degrees counter-clockwise relative to the SVG canvas coordinate system, producing an apparent rotation of the object.

In our module, where we're not going to introduce the concept of nested coordinate systems, let's stick to the first convention in that page: a positive angle rotates the object counter-clockwise within a fixed coordinate system. A right-handed coordinate system. If we agree on that, it's clear that the product of Affine.rotation(angle) * (1.0, 0.0) has been correct but that the matrix representation and its multiplication operation are both wrong (as @perrygeo refers to above).

I'm not sure whether changing both the matrix and the operation while preserving the product is a breaking API change or not, but I'd be fine with releasing a 2.0 version if that helped downstream users. Rasterio, for one, won't be sensitive to this change. Thoughts on that @perrygeo @mwtoews?

@ToddSmall please accept my apologies for how this issue has gone round! Are you still interested in landing this PR?

@ToddSmall
Copy link
Author

@sgillies If you'd prefer to merge your PR, I'm happy to abandon this PR.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants