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

GRASS GRASS user
41

Dylan Beaudette wrote:

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

[snip]

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

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

Right; that definitely isn't correct.

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

42

On Monday 14 August 2006 18:49, Glynn Clements wrote:

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.

Agreed. I was having some terminology issues :slight_smile: . A simple overview of this
interpretation can be found here:

http://local.wasp.uwa.edu.au/~pbourke/colour/bicubic/

Therefore, I would suggest that we name the new version of (incorporating all
three methods) r.resample - with suggestions on how to keep the original
functionality further down:

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

Yes. lets keep the _current_ r.resample functionality, but lets re-name it to
something like r.resamp.default or r.resamp.g or r.resample.g ... (?)
something to denote that this module is mimicking the internal NN resamplng.

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.

Yet another good point Glynn. To keep things simple lets strive for keeping
the concept of r.resample (its future self) consistant across scales.

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

This can be a tricky definition, and one that is a regular source of trouble:
i.e. grid vs pixel registration, etc....

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

just to follow this above comment up: yes. drop r.bilinear, re-name current
r.resample to (?) and create new r.resample based on Gylnn's code - with some
fixes for the bicubic case, as it appears to be broken.

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

Looks like getting the aggregation to work within a C module will be a bit
difficult. r.neighbors was my initial thought for this type of functionality,
but it might complicate the tool chain as you have said down thread...

Cheers,

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

43

On Monday 14 August 2006 19:06, Glynn Clements wrote:

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

Ok. I will accept this. default output to DCELL from the new r.resample sounds
ok to me.

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

44

On Monday 14 August 2006 18:06, Glynn Clements 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;

Some further refinements:

1. when using your new equations, i get correct results when my input raster
actually extends beyond the current region, and when I move from small res to
large res : 10m to 100m.

2. when I try the same, but with a raster that _does not_ extend beyond the
current region, i get odd results, from what looks like accessing an array
with a bad index (?):

(
printed by code [1]
)

[North East]
[4039050.000000 665550.000000]
[R C] [V U]
[0 7] [0 1071644672]
[R0 R1 R2 R3] maprow_f
[-21 -20 -19 -18] -19.500000
0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000
0.000000
73518217792023544336403454714494196428503373132106328150814719669703774177598362215711000331154672582307817565157691
43262327814427811332838925345097477936939540677015455495010019939540775036113228922258615076583985086816808075264.000000
0.00
0000 0.000000
----------------------------------------
c0 c1 c2 c3 :: c
-0.000000 -459488861200147152102521591965588727678146082075664550942591997935648588609989763848193752069716703639423859782235
571453895488401738208302432834068592371058721292313465968438126246221298439757076807641163442286499067926050504704.000000 -0.
000000 -0.000000 :: -25846248442508277305766839548064365931895717116756130990520799883880233109311924216460898553921564579
71759211275075089428162122259777421701184691635832087205307269263246072464460134994803723633557042981544362861557257084034088
96.000000

As this is getting beyond my level of know-how I am going to have to put this
aside for now... I'll take another look later tonight.

Cheers, and thanks for all the pointers

1. code snippet for printing debug info:

printf("[North East]\n");
printf("[%f %f]\n", north, east);

printf("[R C] [V U]\n");
printf("[%i %i] [%i %i]\n", row, col, v, u );

printf("[R0 R1 R2 R3] maprow_f \n");
printf("[%i %i %i %i] %f\n", maprow0, maprow1, maprow2, maprow3, maprow_f );

printf("%f %f %f %f \n", c00, c01, c02, c03);
printf("%f %f %f %f \n", c10, c11, c12, c13);
printf("%f %f %f %f \n", c20, c21, c22, c23);
printf("%f %f %f %f \n", c30, c31, c32, c33);

printf("----------------------------------------\n");
printf("c0 c1 c2 c3 :: c\n");
printf("%f %f %f %f :: %f\n", c0, c1, c2, c3, c);

printf("\n");

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

45

Dylan Beaudette napisa?(a):

Ok. I will accept this. default output to DCELL from the new r.resample sounds
ok to me.

I second that.

Maciek

46

Dylan Beaudette napisa?(a):

just to follow this above comment up: yes. drop r.bilinear, re-name current
r.resample to (?) and create new r.resample based on Gylnn's code

Glynn opts for keeping the old r.resample untouched, and I agree he has
a point. So I won't insist on replacing it, anymore.

Out of other names already suggested, having considered Glynn's
clarifications about interpolation vs resampling, I find r.interpolate
best, in this case.

- with some
fixes for the bicubic case, as it appears to be broken.

Are you refferring to your observations of "weird" slope derivatives of
dems treated with bicubic, ie.

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

?

If so, how does 'gdawalrp -rc' perform for you? Semething like:

$ r.out.gdal input=dem type=Float32 output=dem.tif

$ gdalwarp -te 589980 4913685 609000 4928040 -tr 1 1 -rc dem.tif
gdal_dem.bic.tif
# you need to control the working window with -te switch explicitely to
# mimic Grass region handling

$ r.in.gdal input=gdal_dem.bic.tif output=gdal_dem.bic
$ r.slope.aspect ...

(I'm asking this becaues I always found 'gdalwarp -rc' output correct,
so the eventuall difference might be a pointer.)

???

Maciek

47

On Tuesday 15 August 2006 14:12, Maciej Sieczka wrote:

Dylan Beaudette napisa?(a):
> just to follow this above comment up: yes. drop r.bilinear, re-name
> current r.resample to (?) and create new r.resample based on Gylnn's code

Glynn opts for keeping the old r.resample untouched, and I agree he has
a point. So I won't insist on replacing it, anymore.

good enough for me.

Out of other names already suggested, having considered Glynn's
clarifications about interpolation vs resampling, I find r.interpolate
best, in this case.

ok. Len Coop also suggested 'r.spatial.rescale' thoughts?

> - with some
> fixes for the bicubic case, as it appears to be broken.

Are you refferring to your observations of "weird" slope derivatives of
dems treated with bicubic, ie.

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

?

I was referring to the fact that the bicubic implementation from sample.c (and
revious by Glynn) is not producing the correct output. This is resulting in
both strange output and strange derivative output.

If so, how does 'gdawalrp -rc' perform for you? Semething like:

$ r.out.gdal input=dem type=Float32 output=dem.tif

$ gdalwarp -te 589980 4913685 609000 4928040 -tr 1 1 -rc dem.tif
gdal_dem.bic.tif
# you need to control the working window with -te switch explicitely to
# mimic Grass region handling

$ r.in.gdal input=gdal_dem.bic.tif output=gdal_dem.bic
$ r.slope.aspect ...

(I'm asking this becaues I always found 'gdalwarp -rc' output correct,
so the eventuall difference might be a pointer.)

???

Maciek

For the time being, I am quite sure that the bicubic interpolation is not
working within the currently modified r.bilinear.

However:

1. when using the new equations [1], i get correct results when my input
raster
actually extends beyond the current region, and when I move from small res to
large res : 10m to 100m.

2. when I try the same, but with a raster that _does not_ extend beyond the
current region, i get odd results, from what looks like accessing an array
with a bad index (?):

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

In summary, it looks like we almost have the naming figured out, but there is
still something not quite right with the bicubic method: both the one from
sample.c and the recently corrected version by Glynn.

Cheers,

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

48

Dylan Beaudette wrote:

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

Looks like getting the aggregation to work within a C module will be a bit
difficult. r.neighbors was my initial thought for this type of functionality,
but it might complicate the tool chain as you have said down thread...

The main reason why I suggest r.neighbors is that it already has a set
of aggregate functions.

Hmm; so does r.series. It's been suggested before to make a library of
aggregate functions which could be used by multiple modules. Such
discussions have tended to get bogged down with issues related to
modules such as r.statistics, which need to compute aggregates over a
large number of samples.

If we have 3 modules which need to compute aggregates over a
relatively small number of samples, it's probably worth making a
library from the code for r.series and/or r.neighbors to use with such
modules. We can leave the more complex case (r.statistics) for later.

But first, I need to figure out what's wrong with the updated bicubic
interpolation equations.

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

49

Dylan Beaudette wrote:

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

[snip]

> Try changing the relevant lines to:

[snip]

Some further refinements:

1. when using your new equations, i get correct results when my input raster
actually extends beyond the current region, and when I move from small res to
large res : 10m to 100m.

Is it correct in that both:

a) the interpolated values are close to those obtained by other
methods, and
b) the slope artifacts present in the original version have gone?

?

2. when I try the same, but with a raster that _does not_ extend beyond the
current region, i get odd results, from what looks like accessing an array
with a bad index (?):

Is there any pattern to these errors?

Can you provide a (detailed) recipe which works with spearfish?

Do you get the same errors with the original equations? If only the
equations have changed, a bad index should affect both versions
equally.

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

50

Glynn Clements wrote:

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

[snip]

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

> 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

Right; that definitely isn't correct.

There were two problems with the above:

1. One The 2*c2*u^2 term should have been 4*c2*u^2.

2. The ordering of the indices within the variable names was swapped
(cIJ vs cJI), resulting in each "block" being flipped about the main
diagonal.

I'm confident that the following are correct:

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

This seems to give sane results, including an absence of "grid lines"
in the slope map.

I haven't encountered any totally-out-of-range values, so there may
still be bugs lurking in the I/O code, but I'm fairly sure that the
cubic interpolation equations are finally correct. Well, at least
sensible; there isn't necessarily one "right" way to do cubic
interpolation (although there are plenty wrong ways to do it).

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

51

On Wednesday 16 August 2006 05:11, Glynn Clements wrote:

Glynn Clements wrote:
> > > > 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.
>
> [snip]
>
> > > 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;
> >
> > 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
>
> Right; that definitely isn't correct.

There were two problems with the above:

1. One The 2*c2*u^2 term should have been 4*c2*u^2.

2. The ordering of the indices within the variable names was swapped
(cIJ vs cJI), resulting in each "block" being flipped about the main
diagonal.

I'm confident that the following are correct:

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

This seems to give sane results, including an absence of "grid lines"
in the slope map.

Great. I changed the input, but note that somehow some crazy characters were
included in the in-line text. The results look good for 2 of the 3 test cases
that i have tried. more on this below...

I haven't encountered any totally-out-of-range values, so there may
still be bugs lurking in the I/O code, but I'm fairly sure that the
cubic interpolation equations are finally correct. Well, at least
sensible; there isn't necessarily one "right" way to do cubic
interpolation (although there are plenty wrong ways to do it).

I have summarized three test cases here:
http://169.237.35.250/~dylan/temp/raster_sampling/

1. input raster extent > current region, resampling from 30m to 10m [ok]
2. input raster extent < current region, resampling from 30m to 10m [not ok]
3. input raster extent == current region, resampling from 30m to 10m [ok]

Note that for example 2, the summary from r.info suggests that there is
something wrong with indexing / overflowing :
!!!
min = -0.000000 max = 1356171748364289065433195553253

possibly something to do with dealing with NULL values, and extending the
region???

Cheers,

Dylan

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

52

On Wednesday 16 August 2006 01:20, Glynn Clements wrote:

Dylan Beaudette wrote:
> > 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).
>
> Looks like getting the aggregation to work within a C module will be a
> bit difficult. r.neighbors was my initial thought for this type of
> functionality, but it might complicate the tool chain as you have said
> down thread...

The main reason why I suggest r.neighbors is that it already has a set
of aggregate functions.

Hmm; so does r.series. It's been suggested before to make a library of
aggregate functions which could be used by multiple modules. Such
discussions have tended to get bogged down with issues related to
modules such as r.statistics, which need to compute aggregates over a
large number of samples.

If we have 3 modules which need to compute aggregates over a
relatively small number of samples, it's probably worth making a
library from the code for r.series and/or r.neighbors to use with such
modules. We can leave the more complex case (r.statistics) for later.

Agreed. Lets tackle this after r.interpolate (?) is ready.

But first, I need to figure out what's wrong with the updated bicubic
interpolation equations.

Looks like the bicubic interpolation is almost there
--
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341

53

On Wednesday 16 August 2006 01:26, Glynn Clements wrote:

Dylan Beaudette wrote:
> > > 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.

[snip]

> > Try changing the relevant lines to:

[snip]

> Some further refinements:
>
> 1. when using your new equations, i get correct results when my input
> raster actually extends beyond the current region, and when I move from
> small res to large res : 10m to 100m.

Is it correct in that both:

a) the interpolated values are close to those obtained by other
methods, and
b) the slope artifacts present in the original version have gone?

?

> 2. when I try the same, but with a raster that _does not_ extend beyond
> the current region, i get odd results, from what looks like accessing an
> array with a bad index (?):

Is there any pattern to these errors?

Can you provide a (detailed) recipe which works with spearfish?

Do you get the same errors with the original equations? If only the
equations have changed, a bad index should affect both versions
equally.

I get the same bad results with the original equations as with the new
(correct) equations.

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

54

Dylan Beaudette wrote:

I have summarized three test cases here:
http://169.237.35.250/~dylan/temp/raster_sampling/

1. input raster extent > current region, resampling from 30m to 10m [ok]
2. input raster extent < current region, resampling from 30m to 10m [not ok]
3. input raster extent == current region, resampling from 30m to 10m [ok]

Note that for example 2, the summary from r.info suggests that there is
something wrong with indexing / overflowing :
!!!
min = -0.000000 max = 1356171748364289065433195553253

possibly something to do with dealing with NULL values, and extending the
region???

Yes. I was only considering the bounds of the source map. I've updated
it to set the source region to cover the destination region, plus
neighbours.

I've attached an updated version.

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

(attachments)

main.c.gz (2.05 KB)

55

On Wednesday 16 August 2006 17:28, Glynn Clements wrote:

Dylan Beaudette wrote:
> I have summarized three test cases here:
> http://169.237.35.250/~dylan/temp/raster_sampling/
>
> 1. input raster extent > current region, resampling from 30m to 10m [ok]
> 2. input raster extent < current region, resampling from 30m to 10m [not
> ok] 3. input raster extent == current region, resampling from 30m to 10m
> [ok]
>
> Note that for example 2, the summary from r.info suggests that there is
> something wrong with indexing / overflowing :
> !!!
> min = -0.000000 max = 1356171748364289065433195553253
>
>
> possibly something to do with dealing with NULL values, and extending the
> region???

Yes. I was only considering the bounds of the source map. I've updated
it to set the source region to cover the destination region, plus
neighbours.

I've attached an updated version.

Great! I can confirm that it works as expected. Thanks Glynn for all of the
hard work. Now we just have to decide what to call it, considering previous
conversation of course.

cheers,

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

56

On Wed, 16 Aug 2006, Dylan Beaudette wrote:

Great! I can confirm that it works as expected. Thanks Glynn for all of the
hard work. Now we just have to decide what to call it, considering previous
conversation of course.

To throw in another suggestion, how about r.resamp.interp? This it would mean it would be adjacent in the manuals to the other resampling modules (r.resample and r.resamp.rst) for anyone looking for them, and might also imply that r.resample (the only one with the full word in its name) was kind of a default, basic resampler whilst r.resamp.interp and r.resamp.rst used more advanced and/or interesting resampling methods.

Paul

57

On 8/17/06, Paul Kelly <paul-grass@stjohnspoint.co.uk> wrote:

On Wed, 16 Aug 2006, Dylan Beaudette wrote:

> Great! I can confirm that it works as expected. Thanks Glynn for all of the
> hard work. Now we just have to decide what to call it, considering previous
> conversation of course.

To throw in another suggestion, how about r.resamp.interp? This it would
mean it would be adjacent in the manuals to the other resampling modules
(r.resample and r.resamp.rst) for anyone looking for them, and might also
imply that r.resample (the only one with the full word in its name) was
kind of a default, basic resampler whilst r.resamp.interp and r.resamp.rst
used more advanced and/or interesting resampling methods.

Paul

That sounds good to me.

Cheers,

Dylan

58

On Thu, 17 Aug 2006, Dylan Beaudette wrote:

On 8/17/06, Paul Kelly <paul-grass@stjohnspoint.co.uk> wrote:

To throw in another suggestion, how about r.resamp.interp? This it would
mean it would be adjacent in the manuals to the other resampling modules
(r.resample and r.resamp.rst) for anyone looking for them, and might also
imply that r.resample (the only one with the full word in its name) was
kind of a default, basic resampler whilst r.resamp.interp and r.resamp.rst
used more advanced and/or interesting resampling methods.

That sounds good to me.

Well I would commit it to CVS if it's OK. I briefly edited the r.bilinear man page to create a rough manpage for r.resamp.interp (see below): are there any glaring errors?

Paul

r.resamp.interp(1) Grass User's Manual r.resamp.interp(1)

NAME
        r.resamp.interp - Resamples raster map layers using
        interpolation.

SYNOPSIS
        r.resamp.interp

        r.resamp.interp help

        r.resamp.interp input=name output=name [method=string]
        [--overwrite]

    Flags:
        --overwrite

    Parameters:
        input=name
            Name of input raster map

        output=name
            Name for output raster map

        method=string
            Interpolation method Options: nearest,bilinear,bicubic
            Default: bilinear

DESCRIPTION
        r.resamp.interp fills a grid cell (raster) matrix with
        interpolated values generated from a set of input layer
        data points. A choice of three interpolation methods is
        available; each uses the weighted values of a different
        number of adjacent cells in the input map to determine the
        value of each cell in the output map as follows:

                      nearest neighbour (1 cell)

                      bilinear (4 cells)

                      bicubic (12 cells)

        If there is a current working mask, it applies to the out­
        put raster file. Only those cells falling within the mask
        will be assigned interpolated values. The procedure for
        selection of input data will consider all input data rele­
        vant to interpolating values at the cell centers of the
        current geographic region, ignoring the curent mask. Note
        that for bilinear and bicubic interpolation, cells of the
        output raster that cannot be bounded by the appropriate
        number of input cell centers are set to null.

NOTES
        For longitude-latitude databases, the interpolation algo­
        rithm is based on degree fractions, not on the absolute
        distances between cell centers. Any attempt to implement
        the latter would violate the integrity of the interpola­
        tion method.

SEE ALSO
        g.region, r.resample, r.resamp.rst

AUTHOR
        Glynn Clements

        Last changed: $Date: 2003/04/30 15:17:25 $

        Full index

GRASS 6.3.cvs r.resamp.interp(1)

59

On Thursday 17 August 2006 09:37, Paul Kelly wrote:

On Thu, 17 Aug 2006, Dylan Beaudette wrote:
> On 8/17/06, Paul Kelly <paul-grass@stjohnspoint.co.uk> wrote:
>> To throw in another suggestion, how about r.resamp.interp? This it would
>> mean it would be adjacent in the manuals to the other resampling modules
>> (r.resample and r.resamp.rst) for anyone looking for them, and might
>> also imply that r.resample (the only one with the full word in its name)
>> was kind of a default, basic resampler whilst r.resamp.interp and
>> r.resamp.rst used more advanced and/or interesting resampling methods.
>
> That sounds good to me.

Well I would commit it to CVS if it's OK. I briefly edited the r.bilinear
man page to create a rough manpage for r.resamp.interp (see below): are
there any glaring errors?

Paul

r.resamp.interp(1) Grass User's Manual r.resamp.interp(1)

NAME
        r.resamp.interp - Resamples raster map layers using
        interpolation.

SYNOPSIS
        r.resamp.interp

        r.resamp.interp help

        r.resamp.interp input=name output=name [method=string]
        [--overwrite]

    Flags:
        --overwrite

    Parameters:
        input=name
            Name of input raster map

        output=name
            Name for output raster map

        method=string
            Interpolation method Options: nearest,bilinear,bicubic
            Default: bilinear

DESCRIPTION
        r.resamp.interp fills a grid cell (raster) matrix with
        interpolated values generated from a set of input layer
        data points. A choice of three interpolation methods is
        available; each uses the weighted values of a different
        number of adjacent cells in the input map to determine the
        value of each cell in the output map as follows:

                      nearest neighbour (1 cell)

                      bilinear (4 cells)

                      bicubic (12 cells)

        If there is a current working mask, it applies to the out­
        put raster file. Only those cells falling within the mask
        will be assigned interpolated values. The procedure for
        selection of input data will consider all input data rele­
        vant to interpolating values at the cell centers of the
        current geographic region, ignoring the curent mask. Note
        that for bilinear and bicubic interpolation, cells of the
        output raster that cannot be bounded by the appropriate
        number of input cell centers are set to null.

NOTES
        For longitude-latitude databases, the interpolation algo­
        rithm is based on degree fractions, not on the absolute
        distances between cell centers. Any attempt to implement
        the latter would violate the integrity of the interpola­
        tion method.

SEE ALSO
        g.region, r.resample, r.resamp.rst

AUTHOR
        Glynn Clements

        Last changed: $Date: 2003/04/30 15:17:25 $

        Full index

GRASS 6.3.cvs r.resamp.interp(1)

That looks good for now. Thanks!

Cheers,

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

60

I've been following this thread and am looking forward to using the new
module.

Here's my suggestion on the names.

The new module should inherit the name of r.resample because that is what it
does. Indeed it uses interpolation methods to accomplish this resampling,
but r.interpolate (to me) suggests something that creates a raster surface
from a set of points.

The old module could be renamed to indicate that it resamples using the same
algorithm used in the on-the-fly resampling of g.region in the display. Here
are several suggestions.

r.resample.region
r.dresample
r.gresample
r.region.resample
r.gregion.resample

I suspect that this older module will eventually become unused if the new
module works well. But Glynn makes a good point why we should keep it around
for now.

Michael
__________________________________________
Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics and Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

From: Dylan Beaudette <dylan.beaudette@gmail.com>
Reply-To: <dylan.beaudette@gmail.com>
Date: Wed, 16 Aug 2006 19:03:38 -0700
To: Glynn Clements <glynn@gclements.plus.com>
Cc: <grassuser@grass.itc.it>
Subject: Re: [GRASS-user] integration of sample.c algorithms into r.resample

On Wednesday 16 August 2006 17:28, Glynn Clements wrote:

Dylan Beaudette wrote:

I have summarized three test cases here:
http://169.237.35.250/~dylan/temp/raster_sampling/

1. input raster extent > current region, resampling from 30m to 10m [ok]
2. input raster extent < current region, resampling from 30m to 10m [not
ok] 3. input raster extent == current region, resampling from 30m to 10m
[ok]

Note that for example 2, the summary from r.info suggests that there is
something wrong with indexing / overflowing :
!!!
min = -0.000000 max = 1356171748364289065433195553253

possibly something to do with dealing with NULL values, and extending the
region???

Yes. I was only considering the bounds of the source map. I've updated
it to set the source region to cover the destination region, plus
neighbours.

I've attached an updated version.

Great! I can confirm that it works as expected. Thanks Glynn for all of the
hard work. Now we just have to decide what to call it, considering previous
conversation of course.

cheers,

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