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

[Dylan Beaudette]

On Wednesday 09 August 2006 19:26, Glynn Clements wrote:
> 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.
>

I see. This is starting to get over my head in terms of both C programming and 
the GRASS C API. Are there any _documented_ examples of this sort of 
caching ?

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

OK.

For what I am trying to do, and for the benefit of other modules using an 
interpolator, would it be worthwhile to re-invent what sample.c does, but 
optimized for an entire raster -- or would it be better to somehow re-write 
sample.c to make it more efficient ?

I was hoping for 2 goals: perhaps conflicting (?)

1. a library function which can perform NN, bilinear, cubic interpolation / 
resampling based on an easting/northing pair. this would be useful for both 
sampling at vector points _and_ hopefully re-sampling an entire raster to a 
new resolution

2. a raster function which will perform one of the three interpolations / 
resapling methods on an input raster, saving to an output raster. ideally 
being as fast as possible.

in doing all of this, reduce the amount of code redundancy and bloat, while 
increasing the amount of in-line documentation.

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

helpful...

> 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:

ok. this sounds like a dramatically different approach to what exists in 
r.resample and r.bilinear. Should I just leave these alone, and try and come 
up with an r.cubic ? OR, would it be a good time to replace the cruft of 
r.bilinear with an optimized (?) version of the algorithm (like your 
example), with a cubic method tacked on for good measure, into a new module?

> 	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);
> 		}

Question: will the above code run as-is ? some of these functions 
like 'PUT_DST_CELL()' seem different from the others I have seen in the 
source...?


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

A good test to add to the first section of a new module.. ?


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

Sure. 


Thanks for the feedback. I am interested in getting this moving, but only if 
it will be of benefit to others as well. a hacked version of an old module 
just for a cubic convolution filter is not worth much to me or the community. 

Further pointers would be welcomed, as my C lack of expertise is starting to 
show !   :) 

Cheers,

-- 
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341