Hello Bryce,<br><br>Thanks for your comprehensive email. I agree with many of the points you brought up:<br><ul><li>in RASTER_setRotation, 'rotation' should be 'theta' for clarity</li><li>"rotation about the X axis" should be removed from all documentation, as that is not what we mean</li>
<li>[x|y]skew are synonymous with [x|y]shear</li></ul>In addition, I'm at FOSS4G this week -- any chance you're attending so we can discuss the other details of your email? I'm more interested in doing this "right" than doing it some hacky way (I'm not a math whiz, and my initial efforts may not have been totally correct -- I'm more than happy to change the calculations to be correct).<br>
<br>Thanks,<br>Zwarg<br><br><div class="gmail_quote">On Sat, Sep 10, 2011 at 11:19 AM, Bryce L Nordgren <span dir="ltr"><<a href="mailto:bnordgren@gmail.com">bnordgren@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div class="im">On Sat, Sep 10, 2011 at 4:01 AM, Bborie Park <<a href="mailto:dustymugs@gmail.com">dustymugs@gmail.com</a>> wrote:<br>
><br>
> Hey David,<br>
><br>
</div><div class="im">> I was wondering if you could add the hard limit of 90 degrees for<br>
> RASTER_setRotation. In writing what will become ST_Intersects(raster,<br>
> raster), I was doing some testing on setting skew and what I realized<br>
> is that as the numbers for skew (X and Y) get larger and larger, you<br>
> get infinitely closer to plus/minus 90 degrees but will never there.<br>
> As such, setting a limit saying that the rotation must be less than 90<br>
> degrees sounds reasonable to me.<br>
<br>
</div>This conversation is starting to disturb me a little. :) Of course an<br>
affine transform can rotate to any angle you please (even 5490.2<br>
degrees if you want. It just plain doesn't care.) If we're finding<br>
that it can't, we're doing something wrong. Clearly, because we're<br>
carefully following GDAL's example, we're processing pre-computed<br>
transform coefficients (such as those stored in a geotiff file)<br>
correctly. This is a different animal than calculating the<br>
coefficients correctly ourselves.<br>
<br>
Starting from the beginning: We want to be able to calculate the<br>
parameters for an affine transform which includes the operations:<br>
scaling, translation, rotation, and skew (shearing). Matrices for<br>
these individual operations are found on<br>
<a href="http://en.wikipedia.org/wiki/Transformation_matrix#Examples_in_2D_graphics" target="_blank">http://en.wikipedia.org/wiki/Transformation_matrix#Examples_in_2D_graphics</a><br>
. What ho! We can combine these individual operations willy nilly by<br>
matrix multiplication. But note that these individual matrices are the<br>
only places where individual coefficients represent meaningful<br>
parameters. Once you start the multiplication, the coefficients become<br>
all jumbled up with terms combined in various ways.<br>
<br>
So, using wikipedia plus a little customization, I've labeled the PURE<br>
coefficients (as opposed to our jumbled coefficients) as follows:<br>
<br>
Sx : scale in the x direction<br>
Sy : scale in the y direction<br>
Tx : translation in the x direction<br>
Ty : translation in the y direction<br>
Kx : shearing parallel to x axis<br>
Ky : shearing parallel to y axis<br>
theta : angle of rotation CLOCKWISE around the origin (not around the<br>
x axis, not around the y axis: around the origin; or if you like,<br>
around the invisible Z axis coming out of the screen and poking you in<br>
the eye.)<br>
<br>
The next step is to define an order. I chose: "Scale followed by<br>
rotation followed by shearing (skew)".<br>
<br>
Multiplying these matrices together, as described here<br>
(<a href="http://en.wikipedia.org/wiki/Transformation_matrix#Composing_and_inverting_transformations" target="_blank">http://en.wikipedia.org/wiki/Transformation_matrix#Composing_and_inverting_transformations</a>)<br>
gives a 2x2 matrix, which is not an affine transform yet. Let's label<br>
the coefficients of this matrix as follows<br>
<br>
| O11 O12 |<br>
| O21 O22 |<br>
<br>
(I apologize for the ascii graphics throughout. Gmail is not using a<br>
monospaced font even in "plain text" mode.)<br>
<br>
So what I get for these coefficients (after multiplying in the order<br>
specified) is:<br>
<br>
O11 = Sx * (cos(theta) + Ky sin(theta))<br>
O12 = Sx * (Kx cos(theta) + sin(theta))<br>
O21 = Sy * (-sin(theta) + Ky cos(theta))<br>
O22 = Sy * (-Kx * sin(theta) + cos(theta))<br>
<br>
And the final bit is to add translation by making this into an affine transform:<br>
<br>
| O11 O12 O13 |<br>
| O21 O22 O23 |<br>
| 0 0 1 |<br>
<br>
where:<br>
<br>
O13 = Tx<br>
O23 = Ty<br>
<br>
To be pedantic, this is used as follows:<br>
<br>
| E | | O11 O12 O13 | | i |<br>
| N | = | O21 O22 O23 | | j |<br>
| 1 | | 0 0 1 | | 1 |<br>
<br>
where:<br>
E = easting<br>
N = northing<br>
i = pixel column<br>
j = pixel row<br>
<br>
<br>
The coefficients map onto our jumbled named coefficients as follows :<br>
<br>
ScaleX = O11<br>
SkewX = O12<br>
OffsetX = O13<br>
SkewY = O21<br>
ScaleY = O22<br>
OffsetY = O23<br>
<br>
The important thing to note is that the things we're rather loosely<br>
calling Scale[XY] and Skew[XY] represent all of Scale, Rotation and<br>
Shearing. This dictates that you cannot set these coefficients without<br>
knowing all three.<br>
<br>
Now let's look at RASTER_setRotation...<br>
<br>
xscale = rt_raster_get_x_scale(raster);<br>
yskew = rt_raster_get_y_skew(raster);<br>
psize = sqrt(xscale*xscale + yskew*yskew);<br>
xscale = psize * cos(rotation);<br>
yskew = psize * sin(rotation);<br>
<br>
Clearly, "rotation" == "theta"; "psize" == "Sx"; "xscale" == "O11";<br>
and "yskew" == "O21". In essence, this function looks up the current<br>
values of the jumbled coefficients, attempts to calculate a value for<br>
Sx, and calculates new values for O11 and O12 which do not involve<br>
"Shearing". Furthermore, it purports to solve for Sx without knowing<br>
the previous values for Sx, Sy, Ky, or theta. I'm fairly comfortable<br>
with the claim that this is mathematically impossible regardless of<br>
the particulars.<br>
<br>
In the context of the situation I set up (scale followed by clockwise<br>
rotation followed by shear), the equation for Sx (psize), above is:<br>
<br>
Sx = sqrt(O11^2 + O21^2) ; or<br>
Sx^2 = O11^2 + O21^2<br>
<br>
Let's test that by plopping in the expressions for O11 and O21. I get:<br>
<br>
Sx^2 =? Sx^2 * (cos^2(theta) + 2 Ky cos(theta) sin(theta) + Ky^2<br>
sin^2(theta) ) +<br>
Sy^2 * (sin^2(theta) - 2 Ky cos(theta) sin(theta) +<br>
Ky^2 cos^2(theta))<br>
<br>
I don't see any way for the right-hand side to reduce to Sx^2, other<br>
than the very special case of theta==0 and Ky==0. So at the very<br>
least, we can say that if my math is correct, RASTER_setRotation()<br>
does not implement coefficient calculations which conform to the model<br>
"scale followed by clockwise rotation followed by shear". Given the<br>
way that these pure coefficients get jumbled up, I regard it as<br>
unlikely that it represents any other ordering of these three<br>
operations. So the question is, what does it represent? Where did this<br>
equation in RASTER_setRotation come from?<br>
<br>
Now for the recommendations:<br>
<br>
1] Treat access to the jumbled coefficients of the affine<br>
transformation separately from access to the pure coefficients which<br>
represent something conceptually meaningful. Get rid of the misleading<br>
names for the parameters. Expunge the phrase "rotation about the X<br>
axis" from all documentation unless you REALLY MEAN "rotation in the<br>
Y-Z plane" which you don't. ever. mean.<br>
2] Store the "pure" coefficients with the raster.<br>
3] Explicitly adopt an "orthodox transformation sequence", such as<br>
scale followed by clockwise rotation followed by shearing. If the<br>
geotiff standard specifies an order, go with that. Does GDAL specify<br>
an orthodox model? How does gdal handle the calculation of these<br>
jumbled coefficients?<br>
4] Implement round-trip calculations between the set of "pure"<br>
coefficients and the set of "jumbled" coefficients. The pure->jumbled<br>
calculation is given above (if you like my order of operatons;<br>
otherwise, it's an ELFS). The jumbled->pure calculation may not be<br>
possible, as there are only two equations but there are seven unknowns<br>
(five if you count the fact that you can get Tx and Ty by inspection;<br>
three if you can abandon "shearing").<br>
<br>
w.r.t. 4, you may be able to get two equations per known<br>
geocoordinate. This may be four equations if you know the upper left<br>
and lower right pixel...<br>
<br>
<br>
Food for thought,<br>
<font color="#888888">Bryce<br>
</font><div><div></div><div class="h5">_______________________________________________<br>
postgis-devel mailing list<br>
<a href="mailto:postgis-devel@postgis.refractions.net">postgis-devel@postgis.refractions.net</a><br>
<a href="http://postgis.refractions.net/mailman/listinfo/postgis-devel" target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-devel</a><br>
</div></div></blockquote></div><br>