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

[Glynn Clements]

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 at gclements.plus.com>