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

[Dylan Beaudette]

On 8/14/06, Glynn Clements <glynn at 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 at 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