[postgis-devel] [raster] Implementations of rotation and pixel size.

[Bryce L Nordgren]

On Sat, Sep 10, 2011 at 4:01 AM, Bborie Park <dustymugs at gmail.com> wrote:
>
> Hey David,
>
> I was wondering if you could add the hard limit of 90 degrees for
> RASTER_setRotation.  In writing what will become ST_Intersects(raster,
> raster), I was doing some testing on setting skew and what I realized
> is that as the numbers for skew (X and Y) get larger and larger, you
> get infinitely closer to plus/minus 90 degrees but will never there.
> As such, setting a limit saying that the rotation must be less than 90
> degrees sounds reasonable to me.

This conversation is starting to disturb me a little. :) Of course an
affine transform can rotate to any angle you please (even 5490.2
degrees if you want. It just plain doesn't care.) If we're finding
that it can't, we're doing something wrong. Clearly, because we're
carefully following GDAL's example, we're processing pre-computed
transform coefficients (such as those stored in a geotiff file)
correctly. This is a different animal than calculating the
coefficients correctly ourselves.

Starting from the beginning: We want to be able to calculate the
parameters for an affine transform which includes the operations:
scaling, translation, rotation, and skew (shearing). Matrices for
these individual operations are found on
http://en.wikipedia.org/wiki/Transformation_matrix#Examples_in_2D_graphics
. What ho! We can combine these individual operations willy nilly by
matrix multiplication. But note that these individual matrices are the
only places where individual coefficients represent meaningful
parameters. Once you start the multiplication, the coefficients become
all jumbled up with terms combined in various ways.

So, using wikipedia plus a little customization, I've labeled the PURE
coefficients (as opposed to our jumbled coefficients) as follows:

Sx : scale in the x direction
Sy : scale in the y direction
Tx : translation in the x direction
Ty : translation in the y direction
Kx : shearing parallel to x axis
Ky : shearing parallel to y axis
theta : angle of rotation CLOCKWISE around the origin (not around the
x axis, not around the y axis: around the origin; or if you like,
around the invisible Z axis coming out of the screen and poking you in
the eye.)

The next step is to define an order. I chose: "Scale followed by
rotation followed by shearing (skew)".

Multiplying these matrices together, as described here
(http://en.wikipedia.org/wiki/Transformation_matrix#Composing_and_inverting_transformations)
gives a 2x2 matrix, which is not an affine transform yet. Let's label
the coefficients of this matrix as follows

| O11  O12 |
| O21  O22 |

(I apologize for the ascii graphics throughout. Gmail is not using a
monospaced font even in "plain text" mode.)

So what I get for these coefficients (after multiplying in the order
specified) is:

O11 = Sx * (cos(theta) + Ky sin(theta))
O12 = Sx * (Kx cos(theta) + sin(theta))
O21 = Sy * (-sin(theta) + Ky cos(theta))
O22 = Sy * (-Kx * sin(theta) + cos(theta))

And the final bit is to add translation by making this into an affine transform:

| O11 O12 O13 |
| O21 O22 O23 |
|   0      0     1    |

where:

O13 = Tx
O23 = Ty

To be pedantic, this is used as follows:

| E |         | O11 O12 O13 |  | i |
| N |   =   | O21 O22 O23 |  | j |
| 1 |         |    0     0     1   |   | 1 |

where:
E  = easting
N = northing
i   = pixel column
j   = pixel row


The coefficients map onto our jumbled named coefficients as follows :

ScaleX = O11
SkewX = O12
OffsetX = O13
SkewY = O21
ScaleY = O22
OffsetY = O23

The important thing to note is that the things we're rather loosely
calling Scale[XY] and Skew[XY] represent all of Scale, Rotation and
Shearing. This dictates that you cannot set these coefficients without
knowing all three.

Now let's look at RASTER_setRotation...

    xscale = rt_raster_get_x_scale(raster);
    yskew = rt_raster_get_y_skew(raster);
    psize = sqrt(xscale*xscale + yskew*yskew);
    xscale = psize * cos(rotation);
    yskew = psize * sin(rotation);

Clearly, "rotation" == "theta"; "psize" == "Sx"; "xscale" == "O11";
and "yskew" == "O21". In essence, this function looks up the current
values of the jumbled coefficients, attempts to calculate a value for
Sx, and calculates new values for O11 and O12 which do not involve
"Shearing". Furthermore, it purports to solve for Sx without knowing
the previous values for Sx, Sy, Ky, or theta. I'm fairly comfortable
with the claim that this is mathematically impossible regardless of
the particulars.

In the context of the situation I set up (scale followed by clockwise
rotation followed by shear), the equation for Sx (psize), above is:

Sx = sqrt(O11^2 + O21^2) ; or
Sx^2 = O11^2 + O21^2

Let's test that by plopping in the expressions for O11 and O21. I get:

Sx^2 =? Sx^2 * (cos^2(theta) + 2 Ky cos(theta) sin(theta) + Ky^2
sin^2(theta) ) +
               Sy^2 * (sin^2(theta) - 2 Ky cos(theta) sin(theta) +
Ky^2 cos^2(theta))

I don't see any way for the right-hand side to reduce to Sx^2, other
than the very special case of theta==0 and Ky==0. So at the very
least, we can say that if my math is correct, RASTER_setRotation()
does not implement coefficient calculations which conform to the model
"scale followed by clockwise rotation followed by shear". Given the
way that these pure coefficients get jumbled up, I regard it as
unlikely that it represents any other ordering of these three
operations. So the question is, what does it represent? Where did this
equation in RASTER_setRotation come from?

Now for the recommendations:

1] Treat access to the jumbled coefficients of the affine
transformation separately from access to the pure coefficients which
represent something conceptually meaningful. Get rid of the misleading
names for the parameters. Expunge the phrase "rotation about the X
axis" from all documentation unless you REALLY MEAN "rotation in the
Y-Z plane" which you don't. ever. mean.
2] Store the "pure" coefficients with the raster.
3] Explicitly adopt an "orthodox transformation sequence", such as
scale followed by clockwise rotation followed by shearing. If the
geotiff standard specifies an order, go with that. Does GDAL specify
an orthodox model? How does gdal handle the calculation of these
jumbled coefficients?
4] Implement round-trip calculations between the set of "pure"
coefficients and the set of "jumbled" coefficients. The pure->jumbled
calculation is given above (if you like my order of operatons;
otherwise, it's an ELFS). The jumbled->pure calculation may not be
possible, as there are only two equations but there are seven unknowns
(five if you count the fact that you can get Tx and Ty by inspection;
three if you can abandon "shearing").

w.r.t. 4, you may be able to get two equations per known
geocoordinate. This may be four equations if you know the upper left
and lower right pixel...


Food for thought,
Bryce