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

[Dylan Beaudette]

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