[GRASS-user] integration of sample.c algorithms into r.resample

[Glynn Clements]

Dylan Beaudette wrote:

> > > > If you want any other conversion, you have to keep switching regions,
> > > > so that input maps are read at their native resolution, while output
> > > > maps are written at the desired resolution.
> > >
> > > Yikes. That does complicate things quite a bit. I wonder if
> > > G_get_raster_row() could be modified such that it can perform these three
> > > resampling strategies as well. i.e. such that if one wanted a different
> > > resampling approach in any function that uses G_get_raster_row(), it
> > > would be possible to specify bilinear or cubic convolution... Although
> > > this might be a monumental task.
> >
> > I don't know about "monumental", but it would certainly complicate
> > matters. Apart from anything else, interpolation requires keeping
> > multiple rows in memory.
> >
> > Aside from implementation issues, there are semantic issues, e.g. how
> > you define the interpolation algorithm near the boundaries of the map
> > or if adjoining cells are null.
> >
> > Also, interpolation is only meaningful for scalar quantities, not for
> > discrete categories, and there's no automatic way to distinguish the
> > two (at least for CELL maps; floating-point maps should normally
> > represent scalar quantities, although they could be category maps
> > which have been processed by a module which uses DCELL throughout for
> > convenience).
> 
> Good points. G_get_raster_row() should probably be left alone then. 
> 
> In the mean time, I was hoping to get a general purpose resampling module 
> together mostly for access to a cubic convolution interpolation. 
> 
> In modifying r.bilinear, I have found that:
> 
> 1. the output is not quite correct when the region is not aligned to the input 
> raster, including edge effects
> 2.  the original r.bilinear algorithm is an order of magnitude faster than the 
> algorithm in sample.c

The first thing I notice is that sample.c calls G_get_d_raster_row()
1, 2 or 4 times for each call, which will be a massive performance
hit.

The raster library caches the current row of *raw* data. This
elimnates the read() and decompression, but the rescaling, CELL->DCELL
conversion, embedding nulls etc will be performed on each call. 
Moreover, the raster library only caches a single row, so the bilinear
and bicubic methods will actually perform the read and decompression
steps every time.

Any program which resamples an entire map should be caching the
processed row data, either using the rowio library or doing it itself.

> touching on issue #1 : the output from sample.c , as implemented in this 
> modified version of r.bilinear produces results that are very close to the 
> original r.bilinear when the region is aligned to the input raster. I have 
> attached some additional figures demonstrating this case. 
> 
> Perhaps then, a module based on a generic raster module framework which calls 
> sample.c would do the trick. However, I am not quite sure how to account for 
> edge effects, or a region that is not aligned with the input raster-- the 
> code in r.bilinear is a bit too terse for me to grasp. 
> 
> Any thoughts / insight would be very helpful.

The main point is: don't use sample.c.

Second, read the map at its native resolution, with a slightly
expanded region to allow for the neighbouring cells
(G_get_raster_row() etc will automatically set cells outside the
raster to null). That way, you don't have to write separate cases to
deal with "edge" cells; the handling of null cells will cover it.

Then, it's just a case of iterating over the output grid, translating
the centre of each output cell to a map coordinate which is translated
back to cell coordinates relative to the input map. E.g. for bilinear:

	for (dst_row = 0; dst_row < dst_wind.rows; dst_row++)
		for (dst_col = 0; dst_col < dst_wind.cols; dst_col++)
		{
			double geo_x = dst_wind.west  + (dst_col + 0.5) * dst_wind.ew_res;
			double geo_y = dst_wind.north - (dst_row + 0.5) * dst_wind.ns_res;

			double src_u = (geo_x - src_wind.west ) / src_wind.ew_res - 0.5;
			double src_v = (src_wind.north - geo_y) / src_wind.ns_res - 0.5;

			int col0 = (int) floor(src_u);
			int row0 = (int) floor(src_v);
			int col1 = col0 + 1;
			int row1 = row0 + 1;
			double u_frac = src_u - col0;
			double v_frac = src_v - row0;

			DCELL c00 = GET_SRC_CELL(row0, col0);
			DCELL c01 = GET_SRC_CELL(row0, col1);
			DCELL c10 = GET_SRC_CELL(row1, col0);
			DCELL c11 = GET_SRC_CELL(row1, col1);

			DCELL c0 = c00 * (1 - u_frac) + c01 * u_frac;
			DCELL c1 = c10 * (1 - u_frac) + c11 * u_frac;

			DCELL c = c0 * (1 - v_frac) + c1 * v_frac;

			PUT_DST_CELL(dst_row, dst_col, c);
		}

Note that if src_wind == dst_wind, the dst->geo and geo->src
transformations are exact inverses, and the above reduces to a simple
copy (col0 == dst_col, row0 == dst_row, u_frac == v_frac == 0).

In practice, the code will be slightly complicated by the need to test
for null cells, setting the result to null if any of the four input
cells are null.

-- 
Glynn Clements <glynn at gclements.plus.com>