[GRASS-user] integration of sample.c algorithms into r.resample
glynn at gclements.plus.com
Wed Aug 16 08:11:12 EDT 2006
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.
> > > 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://184.108.40.206/~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
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>
More information about the grass-user