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

[David Zwarg]

Hello Bryce,

Thanks for your comprehensive email. I agree with many of the points you
brought up:

   - in RASTER_setRotation, 'rotation' should be 'theta' for clarity
   - "rotation about the X axis" should be removed from all documentation,
   as that is not what we mean
   - [x|y]skew are synonymous with [x|y]shear

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).

Thanks,
Zwarg

On Sat, Sep 10, 2011 at 11:19 AM, Bryce L Nordgren <bnordgren at gmail.com>wrote:

> 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
> _______________________________________________
> postgis-devel mailing list
> postgis-devel at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-devel/attachments/20110912/5952b3e5/attachment.html>