= How to calculate "real" parameters from affine transformation =

== Introduction == 

This page discusses how to compute physically significant parameters from an arbitrary linear transformation. The reverse process, calculating the affine parameters from the physically significant parameters, is articulated on [wiki:DevWikiAffineParameters another wiki page]. The model used to describe how the raster maps onto real world coordinates is illustrated in the following figure:

[[Image(model.png)]]

This model contains all of the parameters except for translation. Translation may be added in after the coordinates have been scaled, rotated and sheared. The above illustration contains the physically significant parameters which will be calculated by this method. These are: θ,,i,,, θ,,ij,,, and the pixel size along the '''i,,b,,''' and '''j,,b,,''' basis vectors. 

The inputs to this method are the four non-translation (non-offset) parameters of a linear transform.  These are shown below: 

[[Image(geotransform.png)]]

Although these parameters have "common names" within the GIS community, the names are both lengthy and misleading. Within the context of this page, the coefficients o,,11,,, o,,12,,, o,,21,,, and o,,22,, will be used.

This page is divided into sections. Each section describes the computation of one physically significant parameter.

== Pixel size in the '''i''' direction ==

To calculate the pixel size in the '''i''' direction, the '''i''' unit vector is projected using the given transform. This gives the basis vector '''i,,b,,''', the magnitude of which is the pixel size. The basis vector is expressed as follows:

[[Image(basisvector_i.png)]]

The pixel size is then calculated as follows:

[[Image(basisvectormag_i.png)]]

== Pixel size in the '''j''' direction ==

The pixel size in the '''j''' direction is computed in a manner similar to the pixel size in the '''i''' direction. The only difference is that the '''j''' unit vector is used in place of the '''i''' unit vector. This gives the basis vector '''j,,b,,''', the magnitude of which is the pixel size. The basis vector is expressed as follows:

[[Image(basisvector_j.png)]]

The pixel size is then calculated as follows:

[[Image(basisvectormag_j.png)]]

== Rotation == 

The grid is rotated by the angle θ,,i,,. This is the angle between the '''x''' axis of the reference frame and the '''i,,b,,''' basis vector. The angle θ,,i,, is considered positive in the clockwise direction, for consistency with compass headings. The calculation of θ,,i,, is a two-step process: first the magnitude is calculated, then the sign is determined. Both steps involve using the [http://en.wikipedia.org/wiki/Dot_product dot product] to determine the angle between θ,,i,, and one of the axes of the target coordinate system (either the '''x''' axis or the '''y''' axis. The equations in this section refer to the angles and vectors defined in the following figure:

[[Image(calc_theta_i.png)]]

The first step is to calculate the magnitude of the angle between '''i,,b,,''' and the x axis. This is the magnitude of θ,,i,,.

[[Image(thetamag_i.png)]]

The angle θ,,i,, is defined as the angle ''from'' the '''x''' axis ''to'' '''i,,b,,'''. It is positive in the clockwise direction. In the situation described in the above figure, this means that θ,,i,, is negative if '''i,,b,,''' is above the '''x''' axis, and positive if below. We determine whether '''i,,b,,''' is above or below the '''x''' axis by finding the angle between '''i,,b,,''' and the '''y''' axis. If '''i,,b,,''' and the '''y''' axis are separated by less than 90 degrees, '''i,,b,,''' is above the '''x''' axis, otherwise it is below.

[[Image(thetatest_i.png)]]

So, if θ,,test,, is less than 90, θ,,i,, = - abs(θ,,i,,). Otherwise, θ,,i,, = abs(θ,,i,,).

== Basis vector separation angle ==

In this section, the method to calculate θ,,ij,, is presented. This is similar to the method for the calculation of θ,,i,,, but it is accomplished with respect to the rotated reference frame of '''i,,b,,''' and '''i,,bp,,''' instead of the x and y axes. The figure which represents this setup is as follows: 

[[Image(calc_theta_ij.png)]]

The first step is to calculate the magnitude of θ,,ij,,, the angle between '''i,,b,,''' and '''j,,b,,'''.

[[Image(thetamag_ij.png)]]

Next, we need to determine the sign of θ,,ij,, in a manner similar to how the sign for θ,,i,, was determined. The angle θ,,ij,, always represents the angle ''from'' '''i,,b,,''' ''to'' '''j,,b,,''', and is positive counterclockwise for consistency with a right-handed coordinate system. To do this, we first need to calculate '''i,,bp,,''', which is perpendicular to '''i,,b,,''' and forms a [http://en.wikipedia.org/wiki/Cartesian_coordinate_system#Orientation_and_handedness right-hand coordinate system] with '''i,,b,,'''. Observe that '''i,,bp,,''' is '''i,,b,,''' after a 90 degree counterclockwise rotation.

[[Image(basisvector_ip.png)]]

Now we can determine the size of the angle between '''j,,b,,''' and '''i,,bp,,'''. In this situation, any angle less than 90 degrees means that '''j,,b,,''' is on the same side of '''i,,b,,''' as '''i,,bp,,'''. An angle more than 90 degrees means it lies on the opposite side.

[[Image(thetatest_ip.png)]]

So, if θ,,test,, is more than 90 degrees, θ,,ij,, = - abs(θ,,ij,,). Otherwise, θ,,ij,, = abs(θ,,ij,,).


== See also ==

 * Wikipedia article on the [http://en.wikipedia.org/wiki/Dot_product dot product].
 * Wikipedia article on [http://en.wikipedia.org/wiki/Cartesian_coordinate_system coordinate systems].

== The forward calculation == 

On the wiki page describing [wiki:DevWikiAffineParameters the reverse operation], a transform is developed which supports: 

 1. scaling
 1. clockwise rotation around the origin
 1. shearing parallel to the x axis
 1. shearing parallel to the y axis

For the sake of consistency, we will use this same model and attempt to calculate the conceptually meaningful parameters from the transform coefficients. The above model produced the following equations:

 * a,,11,, = o,,11,, = s,,x,, ( (1 + k,,x,, k,,y,,) cosθ + k,,y,, sinθ )
 * a,,12,, = o,,12,, = s,,x,, ( k,,x,, cosθ + sinθ )
 * a,,13,, = t,,x,,
 * a,,21,, = o,,21,, = s,,y,, ( -(1 + k,,x,, k,,y,,) sinθ + k,,y,, cosθ )
 * a,,22,, = o,,22,, = s,,y,, ( - k,,x,, sinθ + cosθ )
 * a,,23,, = t,,y,,

where: 
 * s,,x,, : scale factor in x direction
 * s,,y,, : scale factor in y direction
 * t,,x,, : offset in x direction
 * t,,y,, : offset in y direction
 * θ : angle of rotation clockwise around origin
 * k,,x,, : shearing parallel to x axis
 * k,,y,, : shearing parallel to y axis

== Solving for the meaningful parameters ==

The objective of this page is to find a way to determine values for s,,x,,, s,,y,,, t,,x,,, t,,y,,, θ, k,,x,,, and k,,y,, if all we know is a,,11,,,a,,12,,, ... a,,23,,. Clearly, t,,x,, and t,,y,, are trivial cases, because we can just use the values for a,,13,, and a,,23,, respectively. This leaves us with these four equations: 

 * a,,11,, = o,,11,, = s,,x,, ( (1 + k,,x,, k,,y,,) cosθ + k,,y,, sinθ )
 * a,,12,, = o,,12,, = s,,x,, ( k,,x,, cosθ + sinθ )
 * a,,21,, = o,,21,, = s,,y,, ( -(1 + k,,x,, k,,y,,) sinθ + k,,y,, cosθ )
 * a,,22,, = o,,22,, = s,,y,, ( - k,,x,, sinθ + cosθ )

Unfortunately, these four equations contain five unknowns: s,,x,,, s,,y,,, θ, k,,x,,, and k,,y,,. This represents an ill-posed problem and forces us to simplify. Let's assume that there is no shearing. We are just going to declare upfront that k,,x,,=0 and k,,y,,=0. Remember that we did that because if we're wrong it will screw everything up. This makes the above into: 

 * a,,11,, = o,,11,, = s,,x,, cosθ 
 * a,,12,, = o,,12,, = s,,x,, sinθ 
 * a,,21,, = o,,21,, = - s,,y,, sinθ 
 * a,,22,, = o,,22,, = s,,y,, cosθ 

Now we're on a roll. We can knock out s,,x,, and s,,y,, very easily: 

 * a,,11,,^2^ + a,,12,,^2^ 
  * = s,,x,,^2^ cos^2^θ + s,,x,,^2^ sin^2^θ
  * = s,,x,,^2^ (cos^2^θ + sin^2^θ)
  * = s,,x,,^2^ (1)
 * a,,21,,^2^ + a,,22,,^2^ 
  * = (- s,,y,,)^2^ sin^2^θ + s,,y,,^2^ cos^2^θ
  * = s,,y,,^2^ (sin^2^θ + cos^2^θ)
  * = s,,y,,^2^ (1)

Rewriting as a "final result", this gives: 

 * s,,x,,^2^ = a,,11,,^2^ + a,,12,,^2^
 * s,,y,,^2^ = a,,21,,^2^ + a,,22,,^2^

== Checking against wikipedia ==

Checking this against the wikipedia page on the [http://en.wikipedia.org/wiki/World_file world file], which lists: 
 * "pixel width" = s,,x,, = sqrt(A^2^ + D^2^) = sqrt(a,,11,,^2^ + a,,21,,^2^)
 * "pixel height" = s,,y,, = sqrt(B^2^ + E^2^) = sqrt(a,,12,,^2^ + a,,22,,^2^)

Clearly, this does not match the answer on our page. The three possibilities are: we made a mistake; they made a mistake; or no one made a mistake, but we're using different models... The world file page does not declare it's assumptions about what operations are performed, or which order they are performed in. The world file page, however, also defines it's parameter A (for us: a,,11,,) as "pixel size in x direction". 

'''NOTE: If the order of the model operations changes, in particular if scaling and rotation are swapped, then these pages match wikipedia.''' However, the math on wikipedia is unsourced. Since it has become clear that order is vital, and it is by no means guaranteed that all software performs uses the same model, it is probably vital to go out and determine which software performs the calculations which way before we do anything else.