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

[Glynn Clements]

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>