[raster] C-level index-to-geocoordinate functions

[bnordgren]

Primarily, this patch adds a geocoordinate-to-index function for raster at the C level.

The attached patch changes the internal structure of rt_raster such that it stores forward and reverse affine transformations between pixel index coordinates and geospatial coordinates. The inverse transform is also stored in "struct rt_raster_t" and is automatically recomputed when the accessors are used to change the affine transform parameters. Tidied up the code in rt_api.c to ensure it uses accessors.

The serialized structure is left unchanged, but it is no longer possible to memcpy between struct rt_raster_t and struct rt_raster_serialized_t.

All tests pass. Patch applies against r7789

[bnordgren]

PS: Patch also adds the definition of "struct rt_band_t" to rt_api.h, for external use.

[Bborie Park]

This patch has been superseded by identical functionality added in r7884.

[bnordgren]

With respect, r7884 does two potentially important things differently:

1. It recomputes the matrix inverse of the affine transform on every single call (as opposed to only when the transform changes.) As this is a per-pixel operation, whatever small penalty this may incur could well be multiplied by a factor of several hundred million.
2. There's a strange mix of delegating the operation to GDAL for the cell-to-geopoint direction, but implementing it ourselves for the geopoint-to-cell operation. Worse, neither the comments nor the error messages indicate that this is a matrix inverse operation and there's no mention of where the equation comes from. (Or that "d" is the matrix determinant instead of the scale.) If the commented-out-test-code ever gets deleted, someone could well spend a very long time just trying to figure out what this code is intending to do when investigating some random future problem.

Can we at least add a comment with a link to here: ​http://mathworld.wolfram.com/MatrixInverse.html ?

[Bborie Park]

Replying to bnordgren:

With respect, r7884 does two potentially important things differently:

1. It recomputes the matrix inverse of the affine transform on every single call (as opposed to only when the transform changes.) As this is a per-pixel operation, whatever small penalty this may incur could well be multiplied by a factor of several hundred million.

Only the commented-out validation code generates the matrix inverse. At minimum, the active code only calculates the denominator of the rearranged equation for division by zero. That denominator could be cached…

1. There's a strange mix of delegating the operation to GDAL for the cell-to-geopoint direction, but implementing it ourselves for the geopoint-to-cell operation. Worse, neither the comments nor the error messages indicate that this is a matrix inverse operation and there's no mention of where the equation comes from. (Or that "d" is the matrix determinant instead of the scale.) If the commented-out-test-code ever gets deleted, someone could well spend a very long time just trying to figure out what this code is intending to do when investigating some random future problem.

Neither geopoint_to_cell or cell_to_geopoint use GDAL. "d" is just the denominator and processed prior to everything else just to test for divisiion by zero.

Can we at least add a comment with a link to here: ​http://mathworld.wolfram.com/MatrixInverse.html ?

[bnordgren]

Replying to dustymugs:

Replying to bnordgren:

With respect, r7884 does two potentially important things differently:

Only the commented-out validation code generates the matrix inverse. At minimum, the active code only calculates the denominator of the rearranged equation for division by zero. That denominator could be cached…

Ummmm. Take another look at ​http://mathworld.wolfram.com/MatrixInverse.html

Calculating a disposable matrix inverse is exactly what you're doing. You're just not caching any of the coefficients…and you're mixing in the matrix-vector multiplication at the same time…If you prefer not to believe that you're calculating an inverse matrix, I can go with that too.

My claim here is that if you did cache (or start computing) all six matrix coefficients for the inverse matrix, the forward equations are absolutely no different than the reverse equations. If you'll humor me by letting me count up the number of operations in your forward and reverse transforms:

Forward transform totals:

• 4 multiplies
• 4 additions

Reverse transform totals:

• 6 multiplies
• 2 divides
• 7 subtractions

Room for improvement (if you can reduce your reverse transformation to the same number of operations as the forward one) :

• (8 mult+div - 4 mult) / 8 mult+div = 50%
• (7 sub - 4 add) / 7 sub = 42%

Now divides and multiplies take much longer than adds and subtracts, and you can cut those in half. Big bonus. On top of that, you'll remove nearly half of the adds/subtracts. Still good.

If it was only going to happen once, not such a big deal. But again, it's per pixel. If an operation takes ten minutes with the current code, it'll take five (or less) after the code is reworked to compute/cache the inverse matrix.

Neither geopoint_to_cell or cell_to_geopoint use GDAL. "d" is just the denominator and processed prior to everything else just to test for divisiion by zero.

Whoops. Missed the commented out bits.

[Bborie Park]

Committed using Bryce's suggestions in r7949.