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

GRASS GRASS user
21

On Sunday 13 August 2006 13:52, Glynn Clements wrote:

Maciej Sieczka wrote:
> > For IEEE single precision, FLT_EPSILON is ~1.2e-7, so a value of
> > around 16-17 would have an absolute error of ~2e-6 just due to
> > rounding error. "r.info slope" indicates that the range of the data is
> > between zero and 52.5, so that error seems reasonable.
> >
> >> (A bit OT: what are the Grass equivalents for Gdal's
> >> Int16,UInt16,UInt32,Int32,Float32,Float64,CInt16,CInt32,CFloat32,CFloa
> >>t64? I'd like to put such note into r.out.gdal manual, it always
> >> puzzles me).
> >
> > I'm not sure there are equivalents.
> >
> > Internally, GRASS has three types: CELL (C "int", typically 32-bit
> > signed integer), FCELL (C "float", typically 32-bit IEEE
> > single-precision floating-point) and DCELL (C "double", typically
> > 64-bit IEEE double-precision floating-point).
>
> OK, thanks for clarification. In that case, what are the *most*
> appropriate GDAL datatype equivalents for Grass's CELL, FCELL and DCELL?

AFAICT, the GDAL types would be Int32, Float32 and Float64
respectively.

> >> I was wondering what name should be given to the module. IMO it should
> >> replace r.resample and r.bilinear should be removed. Is it acceptable
> >> for Grass 6.3?
> >
> > I would suggest keeping r.resample. Its output is the result of GRASS'
> > built-in resampling, and is exactly what any raster module would
> > obtain from reading the map with the current region settings. The
> > output from method=nearest in the module which I sent may differ due
> > to rounding error.
>
> If I do some checking and it shows it doesn't differ, would that change
> your attitude (even the backwards-compatibility could be assured by
> setting the default resample method to nearest neighbor (if that's not
> the case yet))?

It would need to produce exactly the same result for all possible
regions, which isn't something which can be verified empirically.

How about this:
1. keep r.resample and update the documentation accordinly.
2. toss r.bilinear, possibly leave a note somewhere in the docs
3. create a new module 'r.resample.3' (?) for the three methods that it would
do? or 'r.interpolate' (?)

Note: the r.interpolate name would make sense for moving from larger to
smaller grid spacing ( 100m -> 10m), however it wouldn't make as much sense
for moving from smaller to larger grid spacing (10m -> 100m). In this case
something like r.aggregate would make more sense. However this potentially
opens a new can of worms- as there are many ways to aggregate cells... which
might call for a new vector / raster module.

Cheers,

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

22

Dylan Beaudette wrote:

How about this:
1. keep r.resample and update the documentation accordinly.
2. toss r.bilinear, possibly leave a note somewhere in the docs

Or replace it with a script which calls "r.interpolate method=bilinear".

3. create a new module 'r.resample.3' (?) for the three methods that it would
do? or 'r.interpolate' (?)

I prefer r.interpolate.

Note: the r.interpolate name would make sense for moving from larger to
smaller grid spacing ( 100m -> 10m), however it wouldn't make as much sense
for moving from smaller to larger grid spacing (10m -> 100m).

Even in that case, it will normally still perform interpolation. The
only case where no interpolation occurs is when the new resolution is
an exact multiple of the old resolution and the alignment is correct.

In this case
something like r.aggregate would make more sense. However this potentially
opens a new can of worms- as there are many ways to aggregate cells... which
might call for a new vector / raster module.

An aggregation tool would be useful in itself. At present, the
only solutions which I know of involve either r.neighbors or
r.average, both of which are a bit awkward.

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

23

Dylan Beaudette wrote:

When I use your new version of r.bilinear method=nearest on a CELL image, the
output is DCELL. While this can be fixed with r.mapcalc it shouldn't be too
hard to add this to the new module (?) .

It's easy enough to change, although I'm not sure whether it's the
right thing to do. It's arguable that the *only* difference between
the methods should be the way in which the result is calculated. Out
of the nine possible combinations of CELL/FCELL/DCELL input and
nearest/bilinear/bicubic methods, at least 8 of them will produce a
floating-point output map, so consistency would favour doing likewise
for the one remaining combination.

Whilst it would be possible to forcibly convert a CELL result to
FCELL/DCELL with r.mapcalc in cases where that was required, it might
not necessarily be obvious that a conversion is necessary until the
first time you actually use method=nearest and end up having to figure
out what happened.

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

24

Glynn Clements napisa?(a):

Dylan Beaudette wrote:

When I use your new version of r.bilinear method=nearest on a CELL image, the
output is DCELL. While this can be fixed with r.mapcalc it shouldn't be too
hard to add this to the new module (?) .

It's easy enough to change, although I'm not sure whether it's the
right thing to do. It's arguable that the *only* difference between
the methods should be the way in which the result is calculated. Out
of the nine possible combinations of CELL/FCELL/DCELL input and
nearest/bilinear/bicubic methods, at least 8 of them will produce a
floating-point output map, so consistency would favour doing likewise
for the one remaining combination.

Whilst it would be possible to forcibly convert a CELL result to
FCELL/DCELL with r.mapcalc in cases where that was required, it might
not necessarily be obvious that a conversion is necessary until the
first time you actually use method=nearest and end up having to figure
out what happened.

Glynn,

That's right. But there is one situation when the output *will* be CELL
for sure - if the input is CELL and method=nearest. Maybe your resample
module could take that into account?

Other combinations are indeed hard to predict, as the most appropriate
output datatype will also depend on the combinations of output region
extent and input/output cell size and shape. So I guess defaulting to
DCELL is very good idea.

Maciek

25

Glynn Clements napisa?(a):

Dylan Beaudette wrote:

How about this:
1. keep r.resample and update the documentation accordinly.
2. toss r.bilinear, possibly leave a note somewhere in the docs

Or replace it with a script which calls "r.interpolate method=bilinear".

3. create a new module 'r.resample.3' (?) for the three methods that it would
do? or 'r.interpolate' (?)

I prefer r.interpolate.

If I may not agree - *inter* suggests this module will be able to fill
null gaps. This is not the case. It will even make the opposite - make
the gaps bigger in the output for method= bilinear and bicubic, if nulls
are present in the input.

This module is for re-sampling and r.resample is the only good name I
can think of currently.

IMHO: drop r.bilinear, replace r.resample. To provide compatibility with
old r.resample (that performed nearest neighbor only) set the default
method=nearest. This will keep things simple and backward compatible.

If it is necessary to mantain r.bilinear compatibilty, make it a script
that calls r.resample method=bilinear.

What do you think?

Maciek

26

Glynn Clements napisa?(a):

Maciej Sieczka wrote:

Internally, GRASS has three types: CELL (C "int", typically 32-bit
signed integer), FCELL (C "float", typically 32-bit IEEE
single-precision floating-point) and DCELL (C "double", typically
64-bit IEEE double-precision floating-point).

OK, thanks for clarification. In that case, what are the *most*
appropriate GDAL datatype equivalents for Grass's CELL, FCELL and DCELL?

AFAICT, the GDAL types would be Int32, Float32 and Float64
respectively.

OK, thanks. I'd like to put it as a hint in the r.out.gdal manual.

Maciek

27

On 8/14/06, Glynn Clements <glynn@gclements.plus.com> wrote:

Dylan Beaudette wrote:

> How about this:
> 1. keep r.resample and update the documentation accordinly.
> 2. toss r.bilinear, possibly leave a note somewhere in the docs

Or replace it with a script which calls "r.interpolate method=bilinear".

> 3. create a new module 'r.resample.3' (?) for the three methods that it would
> do? or 'r.interpolate' (?)

I prefer r.interpolate.

This works for me.

> Note: the r.interpolate name would make sense for moving from larger to
> smaller grid spacing ( 100m -> 10m), however it wouldn't make as much sense
> for moving from smaller to larger grid spacing (10m -> 100m).

Even in that case, it will normally still perform interpolation. The
only case where no interpolation occurs is when the new resolution is
an exact multiple of the old resolution and the alignment is correct.

I am working up a test case for this, cobbling together an
'r.interpolate' via a series of commands to show the differences to
r.resample.

> In this case
> something like r.aggregate would make more sense. However this potentially
> opens a new can of worms- as there are many ways to aggregate cells... which
> might call for a new vector / raster module.

An aggregation tool would be useful in itself. At present, the
only solutions which I know of involve either r.neighbors or
r.average, both of which are a bit awkward.

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

I will post my results a little later today,

Cheers,

Dylan

28

On 8/14/06, Maciej Sieczka <tutey@o2.pl> wrote:

Glynn Clements napisa?(a):
> Dylan Beaudette wrote:
>
>> When I use your new version of r.bilinear method=nearest on a CELL image, the
>> output is DCELL. While this can be fixed with r.mapcalc it shouldn't be too
>> hard to add this to the new module (?) .
>
> It's easy enough to change, although I'm not sure whether it's the
> right thing to do. It's arguable that the *only* difference between
> the methods should be the way in which the result is calculated. Out
> of the nine possible combinations of CELL/FCELL/DCELL input and
> nearest/bilinear/bicubic methods, at least 8 of them will produce a
> floating-point output map, so consistency would favour doing likewise
> for the one remaining combination.

> Whilst it would be possible to forcibly convert a CELL result to
> FCELL/DCELL with r.mapcalc in cases where that was required, it might
> not necessarily be obvious that a conversion is necessary until the
> first time you actually use method=nearest and end up having to figure
> out what happened.

Glynn,

That's right. But there is one situation when the output *will* be CELL
for sure - if the input is CELL and method=nearest. Maybe your resample
module could take that into account?

Other combinations are indeed hard to predict, as the most appropriate
output datatype will also depend on the combinations of output region
extent and input/output cell size and shape. So I guess defaulting to
DCELL is very good idea.

Maciek

I agree with Maciek, in all cases but one defaulting to DCELL seems
like a good idea. However, it might be overly confusing / complicated
when an input map and NN interpolation produce a DCELL map. However, I
found that truncating the DCELL to CELL values worked perfectly well.

Dylan

29

On 8/14/06, Maciej Sieczka <tutey@o2.pl> wrote:

Glynn Clements napisa?(a):
> Dylan Beaudette wrote:
>
>> How about this:
>> 1. keep r.resample and update the documentation accordinly.
>> 2. toss r.bilinear, possibly leave a note somewhere in the docs
>
> Or replace it with a script which calls "r.interpolate method=bilinear".
>
>> 3. create a new module 'r.resample.3' (?) for the three methods that it would
>> do? or 'r.interpolate' (?)
>
> I prefer r.interpolate.

If I may not agree - *inter* suggests this module will be able to fill
null gaps. This is not the case. It will even make the opposite - make
the gaps bigger in the output for method= bilinear and bicubic, if nulls
are present in the input.

This would normally be the case, but the conept of NULL alters the
operational definition of what interpolation is. Either way there is a
potential for misunderstanding. Indeed the methods involved are
traditional _interpolation_ algorithms. As such, they might not be as
useful for moving from small grid sizes to larger grid sized
(aggregating).

This module is for re-sampling and r.resample is the only good name I
can think of currently.

I also like the name r.resample, as it is a bit more general. perhaps
r.resample could be a script that calls r.interpolate / r.aggregate
based on the shift in resolution.

IMHO: drop r.bilinear, replace r.resample. To provide compatibility with
old r.resample (that performed nearest neighbor only) set the default
method=nearest. This will keep things simple and backward compatible.

This also works, but I think that we will need to address the
aggregation problem soon.

If it is necessary to mantain r.bilinear compatibilty, make it a script
that calls r.resample method=bilinear.

What do you think?

Maciek

We are getting closer... !

Cheers,

Dylan

30

On Mon, Aug 14, 2006 at 08:11:56AM -0700, Dylan Beaudette wrote:

On 8/14/06, Glynn Clements <glynn@gclements.plus.com> wrote:
>Dylan Beaudette wrote:

...

>> In this case
>> something like r.aggregate would make more sense. However this
>potentially
>> opens a new can of worms- as there are many ways to aggregate cells...
>which
>> might call for a new vector / raster module.
>
>An aggregation tool would be useful in itself. At present, the
>only solutions which I know of involve either r.neighbors or
>r.average, both of which are a bit awkward.
>
>--
>Glynn Clements <glynn@gclements.plus.com>

You can also use r.statistics for that, but also a bit awkward.
It would be good to have "real" r.aggregate (and maybe drop
merged modules then).

Markus

31

On Monday 14 August 2006 08:11, Dylan Beaudette wrote:

On 8/14/06, Glynn Clements <glynn@gclements.plus.com> wrote:
> Dylan Beaudette wrote:
> > How about this:
> > 1. keep r.resample and update the documentation accordinly.
> > 2. toss r.bilinear, possibly leave a note somewhere in the docs
>
> Or replace it with a script which calls "r.interpolate method=bilinear".
>
> > 3. create a new module 'r.resample.3' (?) for the three methods that it
> > would do? or 'r.interpolate' (?)
>
> I prefer r.interpolate.

This works for me.

> > Note: the r.interpolate name would make sense for moving from larger to
> > smaller grid spacing ( 100m -> 10m), however it wouldn't make as much
> > sense for moving from smaller to larger grid spacing (10m -> 100m).
>
> Even in that case, it will normally still perform interpolation. The
> only case where no interpolation occurs is when the new resolution is
> an exact multiple of the old resolution and the alignment is correct.

I am working up a test case for this, cobbling together an
'r.interpolate' via a series of commands to show the differences to
r.resample.

ooops, this should have said 'cobbling together a 'r.aggregate' module.

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

32

On Monday 14 August 2006 09:36, Dylan Beaudette wrote:

On Monday 14 August 2006 08:11, Dylan Beaudette wrote:
> On 8/14/06, Glynn Clements <glynn@gclements.plus.com> wrote:
> > Dylan Beaudette wrote:
> > > How about this:
> > > 1. keep r.resample and update the documentation accordinly.
> > > 2. toss r.bilinear, possibly leave a note somewhere in the docs
> >
> > Or replace it with a script which calls "r.interpolate
> > method=bilinear".
> >
> > > 3. create a new module 'r.resample.3' (?) for the three methods that
> > > it would do? or 'r.interpolate' (?)
> >
> > I prefer r.interpolate.
>
> This works for me.
>
> > > Note: the r.interpolate name would make sense for moving from larger
> > > to smaller grid spacing ( 100m -> 10m), however it wouldn't make as
> > > much sense for moving from smaller to larger grid spacing (10m ->
> > > 100m).
> >
> > Even in that case, it will normally still perform interpolation. The
> > only case where no interpolation occurs is when the new resolution is
> > an exact multiple of the old resolution and the alignment is correct.
>
> I am working up a test case for this, cobbling together an
> 'r.interpolate' via a series of commands to show the differences to
> r.resample.

ooops, this should have said 'cobbling together a 'r.aggregate' module.

Here is the link to the simple experiment between r.resample and an
aggrigation solution when converting from small to large cell sizes:

http://casoilresource.lawr.ucdavis.edu/drupal/node/275

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

33

On Saturday 12 August 2006 16:11, Glynn Clements wrote:

Maciej Sieczka wrote:
> > I've attached a modified version of r.bilinear which does
> > nearest-neighbour, bilinear and bicubic interpolation. It's only been
> > lightly tested, but the case where the source and destination regions
> > match (exact copy) has been tested. It could use some testing at more
> > extreme scale factors.
> >
> > It would be possible to speed it up a bit, at the expense of
> > simplicity, by caching the horizontal interpolates for each row.
> > However, even in its present form, it's still likely to be a lot
> > faster than using the functions in sample.c.
>
> I have tested and compared your r.bilinear to gdalwarp output. I took
> care to run gdalwarp and r.bilinear with identical settings.
>
> In case of naearest neigbor the output is identical.
>
> In case of biliner and bicubic there is a minor difference between the 2
> programs output, min=-0.000002 and max=0.000002. Exactly the same in
> case of both methods. But this is propably due to the bug that r.in.gdal
> forces FCELL on output http://intevation.de/rt/webrt?serial_num=4897
> while your r.bilinear outputs DCELL. So I guess it is not important.

For IEEE single precision, FLT_EPSILON is ~1.2e-7, so a value of
around 16-17 would have an absolute error of ~2e-6 just due to
rounding error. "r.info slope" indicates that the range of the data is
between zero and 52.5, so that error seems reasonable.

> (A bit OT: what are the Grass equivalents for Gdal's
> Int16,UInt16,UInt32,Int32,Float32,Float64,CInt16,CInt32,CFloat32,CFloat64
>? I'd like to put such note into r.out.gdal manual, it always puzzles me).

I'm not sure there are equivalents.

Internally, GRASS has three types: CELL (C "int", typically 32-bit
signed integer), FCELL (C "float", typically 32-bit IEEE
single-precision floating-point) and DCELL (C "double", typically
64-bit IEEE double-precision floating-point).

> I was wondering what name should be given to the module. IMO it should
> replace r.resample and r.bilinear should be removed. Is it acceptable
> for Grass 6.3?

I would suggest keeping r.resample. Its output is the result of GRASS'
built-in resampling, and is exactly what any raster module would
obtain from reading the map with the current region settings. The
output from method=nearest in the module which I sent may differ due
to rounding error.

After a bit more testing I can see some interesting differences.

using the customized r.bilinear code previously posted:

#set the region
g.region res=10 -ap

#cutout a sample
r.mapcalc "c = elev_meters"

#"resample" 10m DEM data to 1m

# 1. via r.bilinear method=nearest
r.bilinear in=c out=c_n_1m method=nearest
# 2. via r.bilinear method=bilinear
r.bilinear in=c out=c_bc_1m method=bicubic
# 3. via r.bilinear method=bicubic
r.bilinear in=c out=c_b_1m method=bilinear

#visual comparison of slope computed at 1m
r.slope.aspect elev=c_b_1m slope=c_b_1m.slope
r.slope.aspect elev=c_bc_1m slope=c_bc_1m.slope
r.slope.aspect elev=c_n_1m slope=c_n_1m.slope

#reset the colors for each map
r.colors map=c_n_1m.slope color=gyr
r.colors map=c_b_1m.slope color=gyr
r.colors map=c_bc_1m.slope color=gyr

#see images:
##NN
http://169.237.35.250/~dylan/temp/raster_sampling/c_n_1m.slope.png

##Bilinear
http://169.237.35.250/~dylan/temp/raster_sampling/c_b_1m.slope.png

##Bicubic
http://169.237.35.250/~dylan/temp/raster_sampling/c_bc_1m.slope.png

note that the slope computed from the bicubic algorithm looks very similar to
the NN slope map -- all computed via the new r.bilinear code. I was under the
impression that the cubic output would be similar to the bilinear,
but "smoother"...

Thoughts?

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

34

On Monday 14 August 2006 14:08, Dylan Beaudette wrote:

On Saturday 12 August 2006 16:11, Glynn Clements wrote:
> Maciej Sieczka wrote:
> > > I've attached a modified version of r.bilinear which does
> > > nearest-neighbour, bilinear and bicubic interpolation. It's only been
> > > lightly tested, but the case where the source and destination regions
> > > match (exact copy) has been tested. It could use some testing at more
> > > extreme scale factors.
> > >
> > > It would be possible to speed it up a bit, at the expense of
> > > simplicity, by caching the horizontal interpolates for each row.
> > > However, even in its present form, it's still likely to be a lot
> > > faster than using the functions in sample.c.
> >
> > I have tested and compared your r.bilinear to gdalwarp output. I took
> > care to run gdalwarp and r.bilinear with identical settings.
> >
> > In case of naearest neigbor the output is identical.
> >
> > In case of biliner and bicubic there is a minor difference between the
> > 2 programs output, min=-0.000002 and max=0.000002. Exactly the same in
> > case of both methods. But this is propably due to the bug that
> > r.in.gdal forces FCELL on output
> > http://intevation.de/rt/webrt?serial_num=4897 while your r.bilinear
> > outputs DCELL. So I guess it is not important.
>
> For IEEE single precision, FLT_EPSILON is ~1.2e-7, so a value of
> around 16-17 would have an absolute error of ~2e-6 just due to
> rounding error. "r.info slope" indicates that the range of the data is
> between zero and 52.5, so that error seems reasonable.
>
> > (A bit OT: what are the Grass equivalents for Gdal's
> > Int16,UInt16,UInt32,Int32,Float32,Float64,CInt16,CInt32,CFloat32,CFloat
> >64 ? I'd like to put such note into r.out.gdal manual, it always puzzles
> > me).
>
> I'm not sure there are equivalents.
>
> Internally, GRASS has three types: CELL (C "int", typically 32-bit
> signed integer), FCELL (C "float", typically 32-bit IEEE
> single-precision floating-point) and DCELL (C "double", typically
> 64-bit IEEE double-precision floating-point).
>
> > I was wondering what name should be given to the module. IMO it should
> > replace r.resample and r.bilinear should be removed. Is it acceptable
> > for Grass 6.3?
>
> I would suggest keeping r.resample. Its output is the result of GRASS'
> built-in resampling, and is exactly what any raster module would
> obtain from reading the map with the current region settings. The
> output from method=nearest in the module which I sent may differ due
> to rounding error.

After a bit more testing I can see some interesting differences.

using the customized r.bilinear code previously posted:

#set the region
g.region res=10 -ap

#cutout a sample
r.mapcalc "c = elev_meters"

#"resample" 10m DEM data to 1m

# 1. via r.bilinear method=nearest
r.bilinear in=c out=c_n_1m method=nearest
# 2. via r.bilinear method=bilinear
r.bilinear in=c out=c_bc_1m method=bicubic
# 3. via r.bilinear method=bicubic
r.bilinear in=c out=c_b_1m method=bilinear

#visual comparison of slope computed at 1m
r.slope.aspect elev=c_b_1m slope=c_b_1m.slope
r.slope.aspect elev=c_bc_1m slope=c_bc_1m.slope
r.slope.aspect elev=c_n_1m slope=c_n_1m.slope

#reset the colors for each map
r.colors map=c_n_1m.slope color=gyr
r.colors map=c_b_1m.slope color=gyr
r.colors map=c_bc_1m.slope color=gyr

#see images:
##NN
http://169.237.35.250/~dylan/temp/raster_sampling/c_n_1m.slope.png

##Bilinear
http://169.237.35.250/~dylan/temp/raster_sampling/c_b_1m.slope.png

##Bicubic
http://169.237.35.250/~dylan/temp/raster_sampling/c_bc_1m.slope.png

note that the slope computed from the bicubic algorithm looks very similar
to the NN slope map -- all computed via the new r.bilinear code. I was
under the impression that the cubic output would be similar to the
bilinear, but "smoother"...

Thoughts?

note that i left out a

g.region res=1 -a

before performing the r.bilinear commands.

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

35

Dylan Beaudette wrote:

After a bit more testing I can see some interesting differences.

using the customized r.bilinear code previously posted:

#set the region
g.region res=10 -ap

#cutout a sample
r.mapcalc "c = elev_meters"

#"resample" 10m DEM data to 1m

# 1. via r.bilinear method=nearest
r.bilinear in=c out=c_n_1m method=nearest
# 2. via r.bilinear method=bilinear
r.bilinear in=c out=c_bc_1m method=bicubic
# 3. via r.bilinear method=bicubic
r.bilinear in=c out=c_b_1m method=bilinear

#visual comparison of slope computed at 1m
r.slope.aspect elev=c_b_1m slope=c_b_1m.slope
r.slope.aspect elev=c_bc_1m slope=c_bc_1m.slope
r.slope.aspect elev=c_n_1m slope=c_n_1m.slope

#reset the colors for each map
r.colors map=c_n_1m.slope color=gyr
r.colors map=c_b_1m.slope color=gyr
r.colors map=c_bc_1m.slope color=gyr

#see images:
##NN
http://169.237.35.250/~dylan/temp/raster_sampling/c_n_1m.slope.png

##Bilinear
http://169.237.35.250/~dylan/temp/raster_sampling/c_b_1m.slope.png

##Bicubic
http://169.237.35.250/~dylan/temp/raster_sampling/c_bc_1m.slope.png

note that the slope computed from the bicubic algorithm looks very similar to
the NN slope map -- all computed via the new r.bilinear code. I was under the
impression that the cubic output would be similar to the bilinear,
but "smoother"...

The bicubic interpolation formula was taken directly from
lib/gis/sample.c; it wouldn't surprise me if it was wrong.

Hmm.

c = v * (v * (v * (c3 - c2 + c1 - c0) + (c2 - c3 - 2 * c1 + 2 * c0)) + (c2 - c0)) + c1
  = v^3 * (c3 - c2 + c1 - c0) + v^2 * (c2 - c3 - 2*c1 + 2*c0) + v * (c2 - c0) + c1

at v=0 => c1
at v=1 => c2

OK.

dc/dv = 3 * v^2 * (c3 - c2 + c1 - c0) + 2 * v * (c2 - c3 - 2*c1 + 2*c0) + (c2 - c0)
      = v^2 * (3*c3 - 3*c2 + 3*c1 - 3*c0) + v * (2*c2 - 2*c3 - 4*c1 + 4*c0) + (c2 - c0)

at v=0 => c2 - c0
at v=1 => c3 - c1

Not OK.

The actual gradients should be half that (the distance between c2 and
c0 and between c3 and c1 is 2 grid cells). It's doubling the slope of
the boundary tangents.

We want c = k3 * v^3 + k2 * v^2 + k1 * v + k0 such that:

  c(0) = c1
  c(1) = c2
  c'(0) = (c2 - c0)/2
  c'(1) = (c3 - c1)/2

The defintion of c gives:

  c(0) = k0
    => k0 = c1
  c(1) = k3 + k2 + k1 + k0
    => k3 + k2 + k1 + k0 = c2

Differentiating w.r.t. v giv

  c' = 3*k3 * v^2 + 2*k2 * v + k1

  c'(0) = k1
    => k1 = (c2 - c0)/2
  c'(1) = 3*k3 + 2*k2 + k1
    => 3*k3 + 2*k2 + k1 = (c3 - c1)/2

IOW:

  3*k3 + 2*k2 + k1 = (c3 - c1)/2
  k3 + k2 + k1 + k0 = c2
  k1 = (c2 - c0)/2
  k0 = c1

solving gives:
        
  k0 = c1
  k1 = c2/2 - c0/2
  k2 = -c3/2 + 2*c2 - 5*c1/2 + c0
  k3 = c3/2 - 3*c2/2 + 3*c1/2 - c0/2

Try changing the relevant lines to:

  DCELL c0 = (u * (u * (u * (c30 - 3*c20 + 3*c10 - c00) + (-c30 + 2*c20 - 5*c10 + 2*c00)) + (c20 - c00)) + 2*c10)/2;
  DCELL c1 = (u * (u * (u * (c31 - 3*c21 + 3*c11 - c01) + (-c31 + 2*c21 - 5*c11 + 2*c01)) + (c21 - c01)) + 2*c11)/2;
  DCELL c2 = (u * (u * (u * (c32 - 3*c22 + 3*c12 - c02) + (-c32 + 2*c22 - 5*c12 + 2*c02)) + (c22 - c02)) + 2*c12)/2;
  DCELL c3 = (u * (u * (u * (c33 - 3*c23 + 3*c13 - c03) + (-c33 + 2*c23 - 5*c13 + 2*c03)) + (c23 - c03)) + 2*c13)/2;
  DCELL c = (v * (v * (v * (c3 - 3*c2 + 3*c1 - c0 ) + (-c3 + 2*c2 - 5*c1 + 2*c0 )) + (c2 - c0 )) + 2*c1 )/2;

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

36

Maciej Sieczka wrote:

>> How about this:
>> 1. keep r.resample and update the documentation accordinly.
>> 2. toss r.bilinear, possibly leave a note somewhere in the docs
>
> Or replace it with a script which calls "r.interpolate method=bilinear".
>
>> 3. create a new module 'r.resample.3' (?) for the three methods that it would
>> do? or 'r.interpolate' (?)
>
> I prefer r.interpolate.

If I may not agree - *inter* suggests this module will be able to fill
null gaps.

Not really; r.fillnulls is probably the correct name for that.

This is not the case. It will even make the opposite - make
the gaps bigger in the output for method= bilinear and bicubic, if nulls
are present in the input.

This module is for re-sampling and r.resample is the only good name I
can think of currently.

For the bilinear/bicubic methods, the values in the output map aren't
samples, but interpolates of samples.

IMHO: drop r.bilinear, replace r.resample. To provide compatibility with
old r.resample (that performed nearest neighbor only) set the default
method=nearest. This will keep things simple and backward compatible.

"r.<whatever> method=nearest" isn't compatible with r.resample. Even
if you can get the arithmetic such that it exactly matches the current
behaviour of libgis' built-in resampling, there are no guarantees that
this will remain true in future Or even that it will be true for all
combinations of compilation switches (e.g. -ffast-math etc).

IMHO, there needs to be at least one module which does exactly what
r.resample does now, i.e. writes out *exactly* the same data which any
normal raster module will get from G_get_raster_row() etc.

As r.resample already does this, I propose to leave it untouched.

If it is necessary to mantain r.bilinear compatibilty, make it a script
that calls r.resample method=bilinear.

I don't have any problem with either r.bilinear to a script.

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

37

Dylan Beaudette wrote:

Indeed the methods involved are
traditional _interpolation_ algorithms. As such, they might not be as
useful for moving from small grid sizes to larger grid sized
(aggregating).

When used for downscaling, the result for each cell is still an
interpolate of the 2x2 or 4x4 window of source cells surrounding the
centre of the destination cell.

This isn't necessarily wrong, but it's a completely different process
to aggregation.

> This module is for re-sampling and r.resample is the only good name I
> can think of currently.

I also like the name r.resample, as it is a bit more general. perhaps
r.resample could be a script that calls r.interpolate / r.aggregate
based on the shift in resolution.

The existing r.resample needs to remain, even if it is renamed to
allow something else to use that name.

Also, I'm not sure that automatically choosing between interpolation
and aggregation depending upon the scale factor is sensible. While
aggregation is meaningless for scale factors below 1 (and largely
meaningless for scale factors only marginally above 1), interpolation
is meaningful for both reduction and enlargement.

It all depends upon whether your cell values represent samples at
specific points, or aggregates over a rectangular area. Both are
meaningful.

[There's also the case where the scale factor is < 1 in one direction
and > 1 in the other.]

> IMHO: drop r.bilinear, replace r.resample. To provide compatibility with
> old r.resample (that performed nearest neighbor only) set the default
> method=nearest. This will keep things simple and backward compatible.

This also works, but I think that we will need to address the
aggregation problem soon.

Aggregation should probably be done through a modification of
r.neighbors. Rather than using a sliding window at the current
resolution, you would compute the aggregate over the set of source
cells corresponding to a single destination cell.

As with r.bilinear, you would need to repeatedly switch regions to
read the source at its native resolution while writing the result at
the current resolution.

Unlike r.bilinear, where the window size is determined solely by the
method, an aggregation module would have the window size determined by
the scale factor.

This would require some modification to the I/O strategy, specifically
the number of rows kept in memory would vary at run-time, and would
vary from row to row (the number of source rows and columns
corresponding to each destination cell would vary by 1; e.g. a scale
factor 3.5 would result in a destination cell corresponding to 3x3,
3x4, 4x3 or 4x4 source cells).

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

38

Maciej Sieczka wrote:

>> When I use your new version of r.bilinear method=nearest on a CELL image, the
>> output is DCELL. While this can be fixed with r.mapcalc it shouldn't be too
>> hard to add this to the new module (?) .
>
> It's easy enough to change, although I'm not sure whether it's the
> right thing to do. It's arguable that the *only* difference between
> the methods should be the way in which the result is calculated. Out
> of the nine possible combinations of CELL/FCELL/DCELL input and
> nearest/bilinear/bicubic methods, at least 8 of them will produce a
> floating-point output map, so consistency would favour doing likewise
> for the one remaining combination.

> Whilst it would be possible to forcibly convert a CELL result to
> FCELL/DCELL with r.mapcalc in cases where that was required, it might
> not necessarily be obvious that a conversion is necessary until the
> first time you actually use method=nearest and end up having to figure
> out what happened.

That's right. But there is one situation when the output *will* be CELL
for sure - if the input is CELL and method=nearest.

In that case, the output can be converted to CELL without loss of
precision, but the output isn't actually CELL. This matters if you
feed the result to r.mapcalc, where the result of e.g.
"r.mapcalc out = in / 2" will produce different results for integer
and floating point inputs.

Maybe your resample module could take that into account?

As I said, personally I don't think that's a good idea. It's usually
better to let the user handle "special cases" themselves *if* they
want to, rather than handling them automatically then forcing the user
to try to "undo" them.

As the processing "pipeline" gets longer, special cases tend to
combine exponentially, with the end result that the overall process
becomes impossibly hard to document or even understand.

Other combinations are indeed hard to predict, as the most appropriate
output datatype will also depend on the combinations of output region
extent and input/output cell size and shape. So I guess defaulting to
DCELL is very good idea.

Note that the region parameters are all C "double" values, so anything
that involves arithmentic using the interpolation fractions needs to
produce a DCELL result. For the "nearest" method, you could use the
force the output type to match the input type (i.e. CELL in gives CELL
out, FCELL in gives FCELL out, etc).

But if further methods were added (and there is essentially no limit
to the number of interpolation methods available), every method except
nearest would naturally produce DCELL output. IMHO, it's preferable to
simply always produce DCELL output than to create a special case for
one specific method.

In any case, anyone who wants the output to values to exactly match
the input values for the nearest case may as well use the existing
r.resample (which needs to be retained, even if it is renamed to
allow something else to use that name).

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

39

Dylan Beaudette wrote:

> > Whilst it would be possible to forcibly convert a CELL result to
> > FCELL/DCELL with r.mapcalc in cases where that was required, it might
> > not necessarily be obvious that a conversion is necessary until the
> > first time you actually use method=nearest and end up having to figure
> > out what happened.
>
> That's right. But there is one situation when the output *will* be CELL
> for sure - if the input is CELL and method=nearest. Maybe your resample
> module could take that into account?
>
> Other combinations are indeed hard to predict, as the most appropriate
> output datatype will also depend on the combinations of output region
> extent and input/output cell size and shape. So I guess defaulting to
> DCELL is very good idea.

I agree with Maciek, in all cases but one defaulting to DCELL seems
like a good idea. However, it might be overly confusing / complicated
when an input map and NN interpolation produce a DCELL map. However, I
found that truncating the DCELL to CELL values worked perfectly well.

Er, no. What *would* be confusing is having one specific case behave
differently to every other case.

Just because you *can* make method=nearest behave differently to all
of the other methods, that doesn't mean that you should.

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

40

On 8/14/06, Glynn Clements <glynn@gclements.plus.com> wrote:

Dylan Beaudette wrote:

> After a bit more testing I can see some interesting differences.
>
> using the customized r.bilinear code previously posted:
>
> #set the region
> g.region res=10 -ap
>
> #cutout a sample
> r.mapcalc "c = elev_meters"
>
> #"resample" 10m DEM data to 1m
>
> # 1. via r.bilinear method=nearest
> r.bilinear in=c out=c_n_1m method=nearest
> # 2. via r.bilinear method=bilinear
> r.bilinear in=c out=c_bc_1m method=bicubic
> # 3. via r.bilinear method=bicubic
> r.bilinear in=c out=c_b_1m method=bilinear
>
> #visual comparison of slope computed at 1m
> r.slope.aspect elev=c_b_1m slope=c_b_1m.slope
> r.slope.aspect elev=c_bc_1m slope=c_bc_1m.slope
> r.slope.aspect elev=c_n_1m slope=c_n_1m.slope
>
> #reset the colors for each map
> r.colors map=c_n_1m.slope color=gyr
> r.colors map=c_b_1m.slope color=gyr
> r.colors map=c_bc_1m.slope color=gyr
>
> #see images:
> ##NN
> http://169.237.35.250/~dylan/temp/raster_sampling/c_n_1m.slope.png
>
> ##Bilinear
> http://169.237.35.250/~dylan/temp/raster_sampling/c_b_1m.slope.png
>
> ##Bicubic
> http://169.237.35.250/~dylan/temp/raster_sampling/c_bc_1m.slope.png
>
> note that the slope computed from the bicubic algorithm looks very similar to
> the NN slope map -- all computed via the new r.bilinear code. I was under the
> impression that the cubic output would be similar to the bilinear,
> but "smoother"...

The bicubic interpolation formula was taken directly from
lib/gis/sample.c; it wouldn't surprise me if it was wrong.

Hmm.

c = v * (v * (v * (c3 - c2 + c1 - c0) + (c2 - c3 - 2 * c1 + 2 * c0)) + (c2 - c0)) + c1
  = v^3 * (c3 - c2 + c1 - c0) + v^2 * (c2 - c3 - 2*c1 + 2*c0) + v * (c2 - c0) + c1

at v=0 => c1
at v=1 => c2

OK.

dc/dv = 3 * v^2 * (c3 - c2 + c1 - c0) + 2 * v * (c2 - c3 - 2*c1 + 2*c0) + (c2 - c0)
      = v^2 * (3*c3 - 3*c2 + 3*c1 - 3*c0) + v * (2*c2 - 2*c3 - 4*c1 + 4*c0) + (c2 - c0)

at v=0 => c2 - c0
at v=1 => c3 - c1

Not OK.

The actual gradients should be half that (the distance between c2 and
c0 and between c3 and c1 is 2 grid cells). It's doubling the slope of
the boundary tangents.

We want c = k3 * v^3 + k2 * v^2 + k1 * v + k0 such that:

        c(0) = c1
        c(1) = c2
        c'(0) = (c2 - c0)/2
        c'(1) = (c3 - c1)/2

The defintion of c gives:

        c(0) = k0
                => k0 = c1
        c(1) = k3 + k2 + k1 + k0
                => k3 + k2 + k1 + k0 = c2

Differentiating w.r.t. v giv

        c' = 3*k3 * v^2 + 2*k2 * v + k1

        c'(0) = k1
                => k1 = (c2 - c0)/2
        c'(1) = 3*k3 + 2*k2 + k1
                => 3*k3 + 2*k2 + k1 = (c3 - c1)/2

IOW:

        3*k3 + 2*k2 + k1 = (c3 - c1)/2
        k3 + k2 + k1 + k0 = c2
        k1 = (c2 - c0)/2
        k0 = c1

solving gives:

        k0 = c1
        k1 = c2/2 - c0/2
        k2 = -c3/2 + 2*c2 - 5*c1/2 + c0
        k3 = c3/2 - 3*c2/2 + 3*c1/2 - c0/2

Try changing the relevant lines to:

        DCELL c0 = (u * (u * (u * (c30 - 3*c20 + 3*c10 - c00) + (-c30 + 2*c20 - 5*c10 + 2*c00)) + (c20 - c00)) + 2*c10)/2;
        DCELL c1 = (u * (u * (u * (c31 - 3*c21 + 3*c11 - c01) + (-c31 + 2*c21 - 5*c11 + 2*c01)) + (c21 - c01)) + 2*c11)/2;
        DCELL c2 = (u * (u * (u * (c32 - 3*c22 + 3*c12 - c02) + (-c32 + 2*c22 - 5*c12 + 2*c02)) + (c22 - c02)) + 2*c12)/2;
        DCELL c3 = (u * (u * (u * (c33 - 3*c23 + 3*c13 - c03) + (-c33 + 2*c23 - 5*c13 + 2*c03)) + (c23 - c03)) + 2*c13)/2;
        DCELL c = (v * (v * (v * (c3 - 3*c2 + 3*c1 - c0 ) + (-c3 + 2*c2 - 5*c1 + 2*c0 )) + (c2 - c0 )) + 2*c1 )/2;

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

Thanks for the sample code Glynn. I changed your r.bilinear with the
above, but the output is clearly not correct. I do not have time
tonight to dig into the equations, however here is a link to the
output from the above:

http://169.237.35.250/~dylan/temp/raster_sampling/bicubic_v2_broken.png

I wil take a closer look tomorow.

Dylan